]>
Commit | Line | Data |
---|---|---|
200359e8 L |
1 | /* Copyright (C) 2007 Free Software Foundation, Inc. |
2 | ||
3 | This file is part of GCC. | |
4 | ||
5 | GCC is free software; you can redistribute it and/or modify it under | |
6 | the terms of the GNU General Public License as published by the Free | |
7 | Software Foundation; either version 2, or (at your option) any later | |
8 | version. | |
9 | ||
10 | In addition to the permissions in the GNU General Public License, the | |
11 | Free Software Foundation gives you unlimited permission to link the | |
12 | compiled version of this file into combinations with other programs, | |
13 | and to distribute those combinations without any restriction coming | |
14 | from the use of this file. (The General Public License restrictions | |
15 | do apply in other respects; for example, they cover modification of | |
16 | the file, and distribution when not linked into a combine | |
17 | executable.) | |
18 | ||
19 | GCC is distributed in the hope that it will be useful, but WITHOUT ANY | |
20 | WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
21 | FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
22 | for more details. | |
23 | ||
24 | You should have received a copy of the GNU General Public License | |
25 | along with GCC; see the file COPYING. If not, write to the Free | |
26 | Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA | |
27 | 02110-1301, USA. */ | |
28 | ||
29 | /***************************************************************************** | |
30 | * BID64 divide | |
31 | ***************************************************************************** | |
32 | * | |
33 | * Algorithm description: | |
34 | * | |
35 | * if(coefficient_x<coefficient_y) | |
36 | * p = number_digits(coefficient_y) - number_digits(coefficient_x) | |
37 | * A = coefficient_x*10^p | |
38 | * B = coefficient_y | |
39 | * CA= A*10^(15+j), j=0 for A>=B, 1 otherwise | |
40 | * Q = 0 | |
41 | * else | |
42 | * get Q=(int)(coefficient_x/coefficient_y) | |
43 | * (based on double precision divide) | |
44 | * check for exact divide case | |
45 | * Let R = coefficient_x - Q*coefficient_y | |
46 | * Let m=16-number_digits(Q) | |
47 | * CA=R*10^m, Q=Q*10^m | |
48 | * B = coefficient_y | |
49 | * endif | |
50 | * if (CA<2^64) | |
51 | * Q += CA/B (64-bit unsigned divide) | |
52 | * else | |
53 | * get final Q using double precision divide, followed by 3 integer | |
54 | * iterations | |
55 | * if exact result, eliminate trailing zeros | |
56 | * check for underflow | |
57 | * round coefficient to nearest | |
58 | * | |
59 | ****************************************************************************/ | |
60 | ||
61 | #include "bid_internal.h" | |
62 | ||
63 | extern UINT32 __bid_convert_table[5][128][2]; | |
64 | extern SINT8 __bid_factors[][2]; | |
65 | extern UINT8 __bid_packed_10000_zeros[]; | |
66 | ||
67 | ||
68 | #if DECIMAL_CALL_BY_REFERENCE | |
69 | ||
70 | void | |
71 | __bid64_div (UINT64 * pres, UINT64 * px, | |
72 | UINT64 * | |
73 | py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
74 | _EXC_INFO_PARAM) { | |
75 | UINT64 x, y; | |
76 | #else | |
77 | ||
78 | UINT64 | |
79 | __bid64_div (UINT64 x, | |
80 | UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM | |
81 | _EXC_MASKS_PARAM _EXC_INFO_PARAM) { | |
82 | #endif | |
83 | UINT128 CA, CT; | |
84 | UINT64 sign_x, sign_y, coefficient_x, coefficient_y, A, B, QX, PD; | |
85 | UINT64 A2, Q, Q2, B2, B4, B5, R, T, DU, res; | |
86 | UINT64 valid_x, valid_y; | |
87 | SINT64 D; | |
88 | int_double t_scale, tempq, temp_b; | |
89 | int_float tempx, tempy; | |
90 | double da, db, dq, da_h, da_l; | |
91 | int exponent_x = 0, exponent_y = 0, bin_expon_cx; | |
92 | int diff_expon, ed1, ed2, bin_index; | |
93 | int rmode, amount; | |
94 | int nzeros, i, j, k, d5; | |
95 | UINT32 QX32, tdigit[3], digit, digit_h, digit_low; | |
96 | ||
97 | #if DECIMAL_CALL_BY_REFERENCE | |
98 | #if !DECIMAL_GLOBAL_ROUNDING | |
99 | _IDEC_round rnd_mode = *prnd_mode; | |
100 | #endif | |
101 | x = *px; | |
102 | y = *py; | |
103 | #endif | |
104 | ||
105 | valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x); | |
106 | valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y); | |
107 | ||
108 | // unpack arguments, check for NaN or Infinity | |
109 | if (!valid_x) { | |
110 | // x is Inf. or NaN | |
111 | #ifdef SET_STATUS_FLAGS | |
112 | if ((y & SNAN_MASK64) == SNAN_MASK64) // y is sNaN | |
113 | __set_status_flags (pfpsf, INVALID_EXCEPTION); | |
114 | #endif | |
115 | ||
116 | // test if x is NaN | |
117 | if ((x & NAN_MASK64) == NAN_MASK64) { | |
118 | #ifdef SET_STATUS_FLAGS | |
119 | if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN | |
120 | __set_status_flags (pfpsf, INVALID_EXCEPTION); | |
121 | #endif | |
122 | BID_RETURN (x & QUIET_MASK64); | |
123 | } | |
124 | // x is Infinity? | |
125 | if ((x & INFINITY_MASK64) == INFINITY_MASK64) { | |
126 | // check if y is Inf or NaN | |
127 | if ((y & INFINITY_MASK64) == INFINITY_MASK64) { | |
128 | // y==Inf, return NaN | |
129 | #ifdef SET_STATUS_FLAGS | |
130 | if ((y & NAN_MASK64) == INFINITY_MASK64) // Inf/Inf | |
131 | __set_status_flags (pfpsf, INVALID_EXCEPTION); | |
132 | #endif | |
133 | BID_RETURN (NAN_MASK64); | |
134 | } | |
135 | // otherwise return +/-Inf | |
136 | BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64); | |
137 | } | |
138 | // x==0 | |
139 | if (((y & SPECIAL_ENCODING_MASK64) != SPECIAL_ENCODING_MASK64) | |
140 | && !(y << (64 - 53))) { | |
141 | // y==0 , return NaN | |
142 | #ifdef SET_STATUS_FLAGS | |
143 | __set_status_flags (pfpsf, INVALID_EXCEPTION); | |
144 | #endif | |
145 | BID_RETURN (NAN_MASK64); | |
146 | } | |
147 | if (((y & INFINITY_MASK64) != INFINITY_MASK64)) { | |
148 | if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64) | |
149 | exponent_y = ((UINT32) (y >> 51)) & 0x3ff; | |
150 | else | |
151 | exponent_y = ((UINT32) (y >> 53)) & 0x3ff; | |
152 | sign_y = y & 0x8000000000000000ull; | |
153 | ||
154 | exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS; | |
155 | if (exponent_x > DECIMAL_MAX_EXPON_64) | |
156 | exponent_x = DECIMAL_MAX_EXPON_64; | |
157 | else if (exponent_x < 0) | |
158 | exponent_x = 0; | |
159 | BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53)); | |
160 | } | |
161 | ||
162 | } | |
163 | if (!valid_y) { | |
164 | // y is Inf. or NaN | |
165 | ||
166 | // test if y is NaN | |
167 | if ((y & NAN_MASK64) == NAN_MASK64) { | |
168 | #ifdef SET_STATUS_FLAGS | |
169 | if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN | |
170 | __set_status_flags (pfpsf, INVALID_EXCEPTION); | |
171 | #endif | |
172 | BID_RETURN (y & QUIET_MASK64); | |
173 | } | |
174 | // y is Infinity? | |
175 | if ((y & INFINITY_MASK64) == INFINITY_MASK64) { | |
176 | // return +/-0 | |
177 | BID_RETURN (((x ^ y) & 0x8000000000000000ull)); | |
178 | } | |
179 | // y is 0 | |
180 | #ifdef SET_STATUS_FLAGS | |
181 | __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION); | |
182 | #endif | |
183 | BID_RETURN ((sign_x ^ sign_y) | INFINITY_MASK64); | |
184 | } | |
185 | ||
186 | diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS; | |
187 | ||
188 | if (coefficient_x < coefficient_y) { | |
189 | // get number of decimal digits for c_x, c_y | |
190 | ||
191 | //--- get number of bits in the coefficients of x and y --- | |
192 | tempx.d = (float) coefficient_x; | |
193 | tempy.d = (float) coefficient_y; | |
194 | bin_index = (tempy.i - tempx.i) >> 23; | |
195 | ||
196 | A = coefficient_x * __bid_power10_index_binexp[bin_index]; | |
197 | B = coefficient_y; | |
198 | ||
199 | temp_b.d = (double) B; | |
200 | ||
201 | // compare A, B | |
202 | DU = (A - B) >> 63; | |
203 | ed1 = 15 + (int) DU; | |
204 | ed2 = __bid_estimate_decimal_digits[bin_index] + ed1; | |
205 | T = __bid_power10_table_128[ed1].w[0]; | |
206 | __mul_64x64_to_128 (CA, A, T); | |
207 | ||
208 | Q = 0; | |
209 | diff_expon = diff_expon - ed2; | |
210 | ||
211 | // adjust double precision db, to ensure that later A/B - (int)(da/db) > -1 | |
212 | if (coefficient_y < 0x0020000000000000ull) { | |
213 | temp_b.i += 1; | |
214 | db = temp_b.d; | |
215 | } else | |
216 | db = (double) (B + 2 + (B & 1)); | |
217 | ||
218 | } else { | |
219 | // get c_x/c_y | |
220 | ||
221 | // set last bit before conversion to DP | |
222 | A2 = coefficient_x | 1; | |
223 | da = (double) A2; | |
224 | ||
225 | db = (double) coefficient_y; | |
226 | ||
227 | tempq.d = da / db; | |
228 | Q = (UINT64) tempq.d; | |
229 | ||
230 | R = coefficient_x - coefficient_y * Q; | |
231 | ||
232 | // will use to get number of dec. digits of Q | |
233 | bin_expon_cx = (tempq.i >> 52) - 0x3ff; | |
234 | ||
235 | // R<0 ? | |
236 | D = ((SINT64) R) >> 63; | |
237 | Q += D; | |
238 | R += (coefficient_y & D); | |
239 | ||
240 | // exact result ? | |
241 | if (((SINT64) R) <= 0) { | |
242 | // can have R==-1 for coeff_y==1 | |
243 | res = | |
244 | get_BID64 (sign_x ^ sign_y, diff_expon, (Q + R), rnd_mode, | |
245 | pfpsf); | |
246 | BID_RETURN (res); | |
247 | } | |
248 | // get decimal digits of Q | |
249 | DU = __bid_power10_index_binexp[bin_expon_cx] - Q - 1; | |
250 | DU >>= 63; | |
251 | ||
252 | ed2 = 16 - __bid_estimate_decimal_digits[bin_expon_cx] - (int) DU; | |
253 | ||
254 | T = __bid_power10_table_128[ed2].w[0]; | |
255 | __mul_64x64_to_128 (CA, R, T); | |
256 | B = coefficient_y; | |
257 | ||
258 | Q *= __bid_power10_table_128[ed2].w[0]; | |
259 | diff_expon -= ed2; | |
260 | ||
261 | } | |
262 | ||
263 | if (!CA.w[1]) { | |
264 | Q2 = CA.w[0] / B; | |
265 | B2 = B + B; | |
266 | B4 = B2 + B2; | |
267 | R = CA.w[0] - Q2 * B; | |
268 | Q += Q2; | |
269 | } else { | |
270 | ||
271 | // 2^64 | |
272 | t_scale.i = 0x43f0000000000000ull; | |
273 | // convert CA to DP | |
274 | da_h = CA.w[1]; | |
275 | da_l = CA.w[0]; | |
276 | da = da_h * t_scale.d + da_l; | |
277 | ||
278 | // quotient | |
279 | dq = da / db; | |
280 | Q2 = (UINT64) dq; | |
281 | ||
282 | // get w[0] remainder | |
283 | R = CA.w[0] - Q2 * B; | |
284 | ||
285 | // R<0 ? | |
286 | D = ((SINT64) R) >> 63; | |
287 | Q2 += D; | |
288 | R += (B & D); | |
289 | ||
290 | // now R<6*B | |
291 | ||
292 | // quick divide | |
293 | ||
294 | // 4*B | |
295 | B2 = B + B; | |
296 | B4 = B2 + B2; | |
297 | ||
298 | R = R - B4; | |
299 | // R<0 ? | |
300 | D = ((SINT64) R) >> 63; | |
301 | // restore R if negative | |
302 | R += (B4 & D); | |
303 | Q2 += ((~D) & 4); | |
304 | ||
305 | R = R - B2; | |
306 | // R<0 ? | |
307 | D = ((SINT64) R) >> 63; | |
308 | // restore R if negative | |
309 | R += (B2 & D); | |
310 | Q2 += ((~D) & 2); | |
311 | ||
312 | R = R - B; | |
313 | // R<0 ? | |
314 | D = ((SINT64) R) >> 63; | |
315 | // restore R if negative | |
316 | R += (B & D); | |
317 | Q2 += ((~D) & 1); | |
318 | ||
319 | Q += Q2; | |
320 | } | |
321 | ||
322 | #ifdef SET_STATUS_FLAGS | |
323 | if (R) { | |
324 | // set status flags | |
325 | __set_status_flags (pfpsf, INEXACT_EXCEPTION); | |
326 | } | |
327 | #ifndef LEAVE_TRAILING_ZEROS | |
328 | else | |
329 | #endif | |
330 | #else | |
331 | #ifndef LEAVE_TRAILING_ZEROS | |
332 | if (!R) | |
333 | #endif | |
334 | #endif | |
335 | #ifndef LEAVE_TRAILING_ZEROS | |
336 | { | |
337 | // eliminate trailing zeros | |
338 | ||
339 | // check whether CX, CY are short | |
340 | if ((coefficient_x <= 1024) && (coefficient_y <= 1024)) { | |
341 | i = (int) coefficient_y - 1; | |
342 | j = (int) coefficient_x - 1; | |
343 | // difference in powers of 2 __bid_factors for Y and X | |
344 | nzeros = ed2 - __bid_factors[i][0] + __bid_factors[j][0]; | |
345 | // difference in powers of 5 __bid_factors | |
346 | d5 = ed2 - __bid_factors[i][1] + __bid_factors[j][1]; | |
347 | if (d5 < nzeros) | |
348 | nzeros = d5; | |
349 | ||
350 | __mul_64x64_to_128 (CT, Q, __bid_reciprocals10_64[nzeros]); | |
351 | ||
352 | // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 | |
353 | amount = __bid_short_recip_scale[nzeros]; | |
354 | Q = CT.w[1] >> amount; | |
355 | ||
356 | diff_expon += nzeros; | |
357 | } else { | |
358 | tdigit[0] = Q & 0x3ffffff; | |
359 | tdigit[1] = 0; | |
360 | QX = Q >> 26; | |
361 | QX32 = QX; | |
362 | nzeros = 0; | |
363 | ||
364 | for (j = 0; QX32; j++, QX32 >>= 7) { | |
365 | k = (QX32 & 127); | |
366 | tdigit[0] += __bid_convert_table[j][k][0]; | |
367 | tdigit[1] += __bid_convert_table[j][k][1]; | |
368 | if (tdigit[0] >= 100000000) { | |
369 | tdigit[0] -= 100000000; | |
370 | tdigit[1]++; | |
371 | } | |
372 | } | |
373 | ||
374 | digit = tdigit[0]; | |
375 | if (!digit && !tdigit[1]) | |
376 | nzeros += 16; | |
377 | else { | |
378 | if (!digit) { | |
379 | nzeros += 8; | |
380 | digit = tdigit[1]; | |
381 | } | |
382 | // decompose digit | |
383 | PD = (UINT64) digit *0x068DB8BBull; | |
384 | digit_h = (UINT32) (PD >> 40); | |
385 | digit_low = digit - digit_h * 10000; | |
386 | ||
387 | if (!digit_low) | |
388 | nzeros += 4; | |
389 | else | |
390 | digit_h = digit_low; | |
391 | ||
392 | if (!(digit_h & 1)) | |
393 | nzeros += | |
394 | 3 & (UINT32) (__bid_packed_10000_zeros[digit_h >> 3] >> | |
395 | (digit_h & 7)); | |
396 | } | |
397 | ||
398 | if (nzeros) { | |
399 | __mul_64x64_to_128 (CT, Q, __bid_reciprocals10_64[nzeros]); | |
400 | ||
401 | // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 | |
402 | amount = __bid_short_recip_scale[nzeros]; | |
403 | Q = CT.w[1] >> amount; | |
404 | } | |
405 | diff_expon += nzeros; | |
406 | ||
407 | } | |
408 | if (diff_expon >= 0) { | |
409 | res = | |
410 | fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, Q, | |
411 | rnd_mode, pfpsf); | |
412 | BID_RETURN (res); | |
413 | } | |
414 | } | |
415 | #endif | |
416 | ||
417 | if (diff_expon >= 0) { | |
418 | #ifdef IEEE_ROUND_NEAREST | |
419 | // round to nearest code | |
420 | // R*10 | |
421 | R += R; | |
422 | R = (R << 2) + R; | |
423 | B5 = B4 + B; | |
424 | ||
425 | // compare 10*R to 5*B | |
426 | R = B5 - R; | |
427 | // correction for (R==0 && (Q&1)) | |
428 | R -= (Q & 1); | |
429 | // R<0 ? | |
430 | D = ((UINT64) R) >> 63; | |
431 | Q += D; | |
432 | #else | |
433 | #ifdef IEEE_ROUND_NEAREST_TIES_AWAY | |
434 | // round to nearest code | |
435 | // R*10 | |
436 | R += R; | |
437 | R = (R << 2) + R; | |
438 | B5 = B4 + B; | |
439 | ||
440 | // compare 10*R to 5*B | |
441 | R = B5 - R; | |
442 | // correction for (R==0 && (Q&1)) | |
443 | R -= (Q & 1); | |
444 | // R<0 ? | |
445 | D = ((UINT64) R) >> 63; | |
446 | Q += D; | |
447 | #else | |
448 | rmode = rnd_mode; | |
449 | if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2) | |
450 | rmode = 3 - rmode; | |
451 | switch (rmode) { | |
452 | case 0: // round to nearest code | |
453 | case ROUNDING_TIES_AWAY: | |
454 | // R*10 | |
455 | R += R; | |
456 | R = (R << 2) + R; | |
457 | B5 = B4 + B; | |
458 | // compare 10*R to 5*B | |
459 | R = B5 - R; | |
460 | // correction for (R==0 && (Q&1)) | |
461 | R -= ((Q | (rmode >> 2)) & 1); | |
462 | // R<0 ? | |
463 | D = ((UINT64) R) >> 63; | |
464 | Q += D; | |
465 | break; | |
466 | case ROUNDING_DOWN: | |
467 | case ROUNDING_TO_ZERO: | |
468 | break; | |
469 | default: // rounding up | |
470 | Q++; | |
471 | break; | |
472 | } | |
473 | #endif | |
474 | #endif | |
475 | ||
476 | res = | |
477 | fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, Q, rnd_mode, | |
478 | pfpsf); | |
479 | BID_RETURN (res); | |
480 | } else { | |
481 | // UF occurs | |
482 | ||
483 | #ifdef SET_STATUS_FLAGS | |
484 | if ((diff_expon + 16 < 0)) { | |
485 | // set status flags | |
486 | __set_status_flags (pfpsf, INEXACT_EXCEPTION); | |
487 | } | |
488 | #endif | |
489 | rmode = rnd_mode; | |
490 | res = | |
491 | get_BID64_UF (sign_x ^ sign_y, diff_expon, Q, R, rmode, pfpsf); | |
492 | BID_RETURN (res); | |
493 | ||
494 | } | |
495 | ||
496 | ||
497 | } |