OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgcc/] [config/] [libbid/] [bid128_div.c] - Blame information for rev 734

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 734 jeremybenn
/* Copyright (C) 2007, 2009  Free Software Foundation, Inc.
2
 
3
This file is part of GCC.
4
 
5
GCC is free software; you can redistribute it and/or modify it under
6
the terms of the GNU General Public License as published by the Free
7
Software Foundation; either version 3, or (at your option) any later
8
version.
9
 
10
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11
WARRANTY; without even the implied warranty of MERCHANTABILITY or
12
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
13
for more details.
14
 
15
Under Section 7 of GPL version 3, you are granted additional
16
permissions described in the GCC Runtime Library Exception, version
17
3.1, as published by the Free Software Foundation.
18
 
19
You should have received a copy of the GNU General Public License and
20
a copy of the GCC Runtime Library Exception along with this program;
21
see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
22
<http://www.gnu.org/licenses/>.  */
23
 
24
#define BID_128RES
25
#include "bid_div_macros.h"
26
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
27
#include <fenv.h>
28
 
29
#define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
30
#endif
31
 
32
extern UINT32 convert_table[5][128][2];
33
extern SINT8 factors[][2];
34
extern UINT8 packed_10000_zeros[];
35
 
36
BID128_FUNCTION_ARG2 (bid128_div, x, y)
37
 
38
     UINT256 CA4, CA4r, P256;
39
     UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
40
     UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
41
       valid_y;
42
     int_float fx, fy, f64;
43
     UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
44
     int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
45
       digits_q, amount;
46
     int nzeros, i, j, k, d5;
47
     unsigned rmode;
48
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
49
     fexcept_t binaryflags = 0;
50
#endif
51
 
52
valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
53
 
54
  // unpack arguments, check for NaN or Infinity
55
if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
56
    // test if x is NaN
57
if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
58
#ifdef SET_STATUS_FLAGS
59
  if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull ||      // sNaN
60
      (y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull)
61
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
62
#endif
63
  res.w[1] = (CX.w[1]) & QUIET_MASK64;
64
  res.w[0] = CX.w[0];
65
  BID_RETURN (res);
66
}
67
    // x is Infinity?
68
if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
69
  // check if y is Inf. 
70
  if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
71
    // return NaN 
72
  {
73
#ifdef SET_STATUS_FLAGS
74
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
75
#endif
76
    res.w[1] = 0x7c00000000000000ull;
77
    res.w[0] = 0;
78
    BID_RETURN (res);
79
  }
80
  // y is NaN?
81
  if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull))
82
    // return NaN 
83
  {
84
    // return +/-Inf
85
    res.w[1] = ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) |
86
      0x7800000000000000ull;
87
    res.w[0] = 0;
88
    BID_RETURN (res);
89
  }
90
}
91
    // x is 0
92
if ((y.w[1] & 0x7800000000000000ull) < 0x7800000000000000ull) {
93
  if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) {
94
#ifdef SET_STATUS_FLAGS
95
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
96
#endif
97
    // x=y=0, return NaN
98
    res.w[1] = 0x7c00000000000000ull;
99
    res.w[0] = 0;
100
    BID_RETURN (res);
101
  }
102
  // return 0
103
  res.w[1] = (x.w[1] ^ y.w[1]) & 0x8000000000000000ull;
104
  exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
105
  if (exponent_x > DECIMAL_MAX_EXPON_128)
106
    exponent_x = DECIMAL_MAX_EXPON_128;
107
  else if (exponent_x < 0)
108
    exponent_x = 0;
109
  res.w[1] |= (((UINT64) exponent_x) << 49);
110
  res.w[0] = 0;
111
  BID_RETURN (res);
112
}
113
}
114
if (!valid_y) {
115
  // y is Inf. or NaN
116
 
117
  // test if y is NaN
118
  if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
119
#ifdef SET_STATUS_FLAGS
120
    if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull)      // sNaN
121
      __set_status_flags (pfpsf, INVALID_EXCEPTION);
122
#endif
123
    res.w[1] = CY.w[1] & QUIET_MASK64;
124
    res.w[0] = CY.w[0];
125
    BID_RETURN (res);
126
  }
127
  // y is Infinity?
128
  if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
129
    // return +/-0
130
    res.w[1] = sign_x ^ sign_y;
131
    res.w[0] = 0;
132
    BID_RETURN (res);
133
  }
134
  // y is 0, return +/-Inf
135
#ifdef SET_STATUS_FLAGS
136
  __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
137
#endif
138
  res.w[1] =
139
    ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
140
  res.w[0] = 0;
141
  BID_RETURN (res);
142
}
143
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
144
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
145
#endif
146
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
147
 
148
if (__unsigned_compare_gt_128 (CY, CX)) {
149
  // CX < CY
150
 
151
  // 2^64
152
  f64.i = 0x5f800000;
153
 
154
  // fx ~ CX,   fy ~ CY
155
  fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
156
  fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
157
  // expon_cy - expon_cx
158
  bin_index = (fy.i - fx.i) >> 23;
159
 
160
  if (CX.w[1]) {
161
    T = power10_index_binexp_128[bin_index].w[0];
162
    __mul_64x128_short (CA, T, CX);
163
  } else {
164
    T128 = power10_index_binexp_128[bin_index];
165
    __mul_64x128_short (CA, CX.w[0], T128);
166
  }
167
 
168
  ed2 = 33;
169
  if (__unsigned_compare_gt_128 (CY, CA))
170
    ed2++;
171
 
172
  T128 = power10_table_128[ed2];
173
  __mul_128x128_to_256 (CA4, CA, T128);
174
 
175
  ed2 += estimate_decimal_digits[bin_index];
176
  CQ.w[0] = CQ.w[1] = 0;
177
  diff_expon = diff_expon - ed2;
178
 
179
} else {
180
  // get CQ = CX/CY
181
  __div_128_by_128 (&CQ, &CR, CX, CY);
182
 
183
  if (!CR.w[1] && !CR.w[0]) {
184
    get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
185
                pfpsf);
186
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
187
    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
188
#endif
189
    BID_RETURN (res);
190
  }
191
  // get number of decimal digits in CQ
192
  // 2^64
193
  f64.i = 0x5f800000;
194
  fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
195
  // binary expon. of CQ
196
  bin_expon = (fx.i - 0x3f800000) >> 23;
197
 
198
  digits_q = estimate_decimal_digits[bin_expon];
199
  TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
200
  TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
201
  if (__unsigned_compare_ge_128 (CQ, TP128))
202
    digits_q++;
203
 
204
  ed2 = 34 - digits_q;
205
  T128.w[0] = power10_table_128[ed2].w[0];
206
  T128.w[1] = power10_table_128[ed2].w[1];
207
  __mul_128x128_to_256 (CA4, CR, T128);
208
  diff_expon = diff_expon - ed2;
209
  __mul_128x128_low (CQ, CQ, T128);
210
 
211
}
212
 
213
__div_256_by_128 (&CQ, &CA4, CY);
214
 
215
#ifdef SET_STATUS_FLAGS
216
if (CA4.w[0] || CA4.w[1]) {
217
  // set status flags
218
  __set_status_flags (pfpsf, INEXACT_EXCEPTION);
219
}
220
#ifndef LEAVE_TRAILING_ZEROS
221
else
222
#endif
223
#else
224
#ifndef LEAVE_TRAILING_ZEROS
225
if (!CA4.w[0] && !CA4.w[1])
226
#endif
227
#endif
228
#ifndef LEAVE_TRAILING_ZEROS
229
  // check whether result is exact
230
{
231
  // check whether CX, CY are short
232
  if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
233
    i = (int) CY.w[0] - 1;
234
    j = (int) CX.w[0] - 1;
235
    // difference in powers of 2 factors for Y and X
236
    nzeros = ed2 - factors[i][0] + factors[j][0];
237
    // difference in powers of 5 factors
238
    d5 = ed2 - factors[i][1] + factors[j][1];
239
    if (d5 < nzeros)
240
      nzeros = d5;
241
    // get P*(2^M[extra_digits])/10^extra_digits
242
    __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
243
 
244
    // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
245
    amount = recip_scale[nzeros];
246
    __shr_128_long (CQ, Qh, amount);
247
 
248
    diff_expon += nzeros;
249
  } else {
250
    // decompose Q as Qh*10^17 + Ql
251
    //T128 = reciprocals10_128[17];
252
    T128.w[0] = 0x44909befeb9fad49ull;
253
    T128.w[1] = 0x000b877aa3236a4bull;
254
    __mul_128x128_to_256 (P256, CQ, T128);
255
    //amount = recip_scale[17];
256
    Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
257
    Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
258
 
259
    if (!Q_low) {
260
      diff_expon += 17;
261
 
262
      tdigit[0] = Q_high & 0x3ffffff;
263
      tdigit[1] = 0;
264
      QX = Q_high >> 26;
265
      QX32 = QX;
266
      nzeros = 0;
267
 
268
      for (j = 0; QX32; j++, QX32 >>= 7) {
269
        k = (QX32 & 127);
270
        tdigit[0] += convert_table[j][k][0];
271
        tdigit[1] += convert_table[j][k][1];
272
        if (tdigit[0] >= 100000000) {
273
          tdigit[0] -= 100000000;
274
          tdigit[1]++;
275
        }
276
      }
277
 
278
      if (tdigit[1] >= 100000000) {
279
        tdigit[1] -= 100000000;
280
        if (tdigit[1] >= 100000000)
281
          tdigit[1] -= 100000000;
282
      }
283
 
284
      digit = tdigit[0];
285
      if (!digit && !tdigit[1])
286
        nzeros += 16;
287
      else {
288
        if (!digit) {
289
          nzeros += 8;
290
          digit = tdigit[1];
291
        }
292
        // decompose digit
293
        PD = (UINT64) digit *0x068DB8BBull;
294
        digit_h = (UINT32) (PD >> 40);
295
        digit_low = digit - digit_h * 10000;
296
 
297
        if (!digit_low)
298
          nzeros += 4;
299
        else
300
          digit_h = digit_low;
301
 
302
        if (!(digit_h & 1))
303
          nzeros +=
304
            3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
305
                          (digit_h & 7));
306
      }
307
 
308
      if (nzeros) {
309
        __mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
310
 
311
        // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
312
        amount = short_recip_scale[nzeros];
313
        CQ.w[0] = CQ.w[1] >> amount;
314
      } else
315
        CQ.w[0] = Q_high;
316
      CQ.w[1] = 0;
317
 
318
      diff_expon += nzeros;
319
    } else {
320
      tdigit[0] = Q_low & 0x3ffffff;
321
      tdigit[1] = 0;
322
      QX = Q_low >> 26;
323
      QX32 = QX;
324
      nzeros = 0;
325
 
326
      for (j = 0; QX32; j++, QX32 >>= 7) {
327
        k = (QX32 & 127);
328
        tdigit[0] += convert_table[j][k][0];
329
        tdigit[1] += convert_table[j][k][1];
330
        if (tdigit[0] >= 100000000) {
331
          tdigit[0] -= 100000000;
332
          tdigit[1]++;
333
        }
334
      }
335
 
336
      if (tdigit[1] >= 100000000) {
337
        tdigit[1] -= 100000000;
338
        if (tdigit[1] >= 100000000)
339
          tdigit[1] -= 100000000;
340
      }
341
 
342
      digit = tdigit[0];
343
      if (!digit && !tdigit[1])
344
        nzeros += 16;
345
      else {
346
        if (!digit) {
347
          nzeros += 8;
348
          digit = tdigit[1];
349
        }
350
        // decompose digit
351
        PD = (UINT64) digit *0x068DB8BBull;
352
        digit_h = (UINT32) (PD >> 40);
353
        digit_low = digit - digit_h * 10000;
354
 
355
        if (!digit_low)
356
          nzeros += 4;
357
        else
358
          digit_h = digit_low;
359
 
360
        if (!(digit_h & 1))
361
          nzeros +=
362
            3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
363
                          (digit_h & 7));
364
      }
365
 
366
      if (nzeros) {
367
        // get P*(2^M[extra_digits])/10^extra_digits
368
        __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
369
 
370
        //now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
371
        amount = recip_scale[nzeros];
372
        __shr_128 (CQ, Qh, amount);
373
      }
374
      diff_expon += nzeros;
375
 
376
    }
377
  }
378
  get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
379
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
380
  (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
381
#endif
382
  BID_RETURN (res);
383
}
384
#endif
385
 
386
if (diff_expon >= 0) {
387
#ifdef IEEE_ROUND_NEAREST
388
  // rounding
389
  // 2*CA4 - CY
390
  CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
391
  CA4r.w[0] = CA4.w[0] + CA4.w[0];
392
  __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
393
  CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
394
 
395
  D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
396
  carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
397
 
398
  CQ.w[0] += carry64;
399
  if (CQ.w[0] < carry64)
400
    CQ.w[1]++;
401
#else
402
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
403
  // rounding
404
  // 2*CA4 - CY
405
  CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
406
  CA4r.w[0] = CA4.w[0] + CA4.w[0];
407
  __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
408
  CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
409
 
410
  D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
411
  carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
412
 
413
  CQ.w[0] += carry64;
414
  if (CQ.w[0] < carry64)
415
    CQ.w[1]++;
416
#else
417
  rmode = rnd_mode;
418
  if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
419
    rmode = 3 - rmode;
420
  switch (rmode) {
421
  case ROUNDING_TO_NEAREST:     // round to nearest code
422
    // rounding
423
    // 2*CA4 - CY
424
    CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
425
    CA4r.w[0] = CA4.w[0] + CA4.w[0];
426
    __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
427
    CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
428
    D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
429
    carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
430
    CQ.w[0] += carry64;
431
    if (CQ.w[0] < carry64)
432
      CQ.w[1]++;
433
    break;
434
  case ROUNDING_TIES_AWAY:
435
    // rounding
436
    // 2*CA4 - CY
437
    CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
438
    CA4r.w[0] = CA4.w[0] + CA4.w[0];
439
    __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
440
    CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
441
    D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
442
    carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
443
    CQ.w[0] += carry64;
444
    if (CQ.w[0] < carry64)
445
      CQ.w[1]++;
446
    break;
447
  case ROUNDING_DOWN:
448
  case ROUNDING_TO_ZERO:
449
    break;
450
  default:      // rounding up
451
    CQ.w[0]++;
452
    if (!CQ.w[0])
453
      CQ.w[1]++;
454
    break;
455
  }
456
#endif
457
#endif
458
 
459
} else {
460
#ifdef SET_STATUS_FLAGS
461
  if (CA4.w[0] || CA4.w[1]) {
462
    // set status flags
463
    __set_status_flags (pfpsf, INEXACT_EXCEPTION);
464
  }
465
#endif
466
 
467
  handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
468
                     CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
469
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
470
  (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
471
#endif
472
  BID_RETURN (res);
473
 
474
}
475
 
476
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
477
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
478
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
479
#endif
480
BID_RETURN (res);
481
}
482
 
483
 
484
//#define LEAVE_TRAILING_ZEROS
485
 
486
TYPE0_FUNCTION_ARGTYPE1_ARGTYPE2 (UINT128, bid128dd_div, UINT64, x,
487
                                  UINT64, y)
488
 
489
     UINT256 CA4, CA4r, P256;
490
     UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
491
     UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
492
       valid_y;
493
     int_float fx, fy, f64;
494
     UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
495
     int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
496
       digits_q, amount;
497
     int nzeros, i, j, k, d5;
498
     unsigned rmode;
499
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
500
     fexcept_t binaryflags = 0;
501
#endif
502
 
503
valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], y);
504
 
505
        // unpack arguments, check for NaN or Infinity
506
CX.w[1] = 0;
507
if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], (x))) {
508
#ifdef SET_STATUS_FLAGS
509
if ((y & SNAN_MASK64) == SNAN_MASK64)   // y is sNaN
510
  __set_status_flags (pfpsf, INVALID_EXCEPTION);
511
#endif
512
 
513
    // test if x is NaN
514
if ((x & NAN_MASK64) == NAN_MASK64) {
515
#ifdef SET_STATUS_FLAGS
516
  if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
517
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
518
#endif
519
  res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
520
  __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
521
  res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
522
  BID_RETURN (res);
523
}
524
           // x is Infinity?
525
if (((x) & 0x7800000000000000ull) == 0x7800000000000000ull) {
526
  // check if y is Inf.
527
  if ((((y) & 0x7c00000000000000ull) == 0x7800000000000000ull))
528
    // return NaN 
529
  {
530
#ifdef SET_STATUS_FLAGS
531
  __set_status_flags (pfpsf, INVALID_EXCEPTION);
532
#endif
533
  res.w[1] = 0x7c00000000000000ull;
534
  res.w[0] = 0;
535
    BID_RETURN (res);
536
  }
537
  if ((((y) & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
538
  // otherwise return +/-Inf
539
  res.w[1] =
540
    (((x) ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
541
  res.w[0] = 0;
542
  BID_RETURN (res);
543
  }
544
}
545
           // x is 0
546
if ((((y) & 0x7800000000000000ull) != 0x7800000000000000ull)) {
547
    if(!CY.w[0]) {
548
#ifdef SET_STATUS_FLAGS
549
  __set_status_flags (pfpsf, INVALID_EXCEPTION);
550
#endif
551
  // x=y=0, return NaN
552
  res.w[1] = 0x7c00000000000000ull;
553
  res.w[0] = 0;
554
  BID_RETURN (res);
555
}
556
           // return 0
557
res.w[1] = ((x) ^ (y)) & 0x8000000000000000ull;
558
if (((y) & 0x6000000000000000ull) == 0x6000000000000000ull)
559
  exponent_y = ((UINT32) ((y) >> 51)) & 0x3ff;
560
else
561
  exponent_y = ((UINT32) ((y) >> 53)) & 0x3ff;
562
exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
563
if (exponent_x > DECIMAL_MAX_EXPON_128)
564
  exponent_x = DECIMAL_MAX_EXPON_128;
565
else if (exponent_x < 0)
566
  exponent_x = 0;
567
res.w[1] |= (((UINT64) exponent_x) << 49);
568
res.w[0] = 0;
569
BID_RETURN (res);
570
}
571
}
572
 
573
CY.w[1] = 0;
574
if (!valid_y) {
575
  // y is Inf. or NaN
576
 
577
  // test if y is NaN
578
  if ((y & NAN_MASK64) == NAN_MASK64) {
579
#ifdef SET_STATUS_FLAGS
580
    if ((y & SNAN_MASK64) == SNAN_MASK64)       // sNaN
581
      __set_status_flags (pfpsf, INVALID_EXCEPTION);
582
#endif
583
  res.w[0] = (CY.w[0] & 0x0003ffffffffffffull);
584
  __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
585
  res.w[1] |= ((CY.w[0]) & 0xfc00000000000000ull);
586
    BID_RETURN (res);
587
  }
588
  // y is Infinity?
589
  if (((y) & 0x7800000000000000ull) == 0x7800000000000000ull) {
590
    // return +/-0
591
    res.w[1] = sign_x ^ sign_y;
592
    res.w[0] = 0;
593
    BID_RETURN (res);
594
  }
595
  // y is 0, return +/-Inf
596
  res.w[1] =
597
    (((x) ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
598
  res.w[0] = 0;
599
#ifdef SET_STATUS_FLAGS
600
  __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
601
#endif
602
  BID_RETURN (res);
603
}
604
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
605
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
606
#endif
607
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
608
 
609
if (__unsigned_compare_gt_128 (CY, CX)) {
610
  // CX < CY
611
 
612
  // 2^64
613
  f64.i = 0x5f800000;
614
 
615
  // fx ~ CX,   fy ~ CY
616
  fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
617
  fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
618
  // expon_cy - expon_cx
619
  bin_index = (fy.i - fx.i) >> 23;
620
 
621
  if (CX.w[1]) {
622
    T = power10_index_binexp_128[bin_index].w[0];
623
    __mul_64x128_short (CA, T, CX);
624
  } else {
625
    T128 = power10_index_binexp_128[bin_index];
626
    __mul_64x128_short (CA, CX.w[0], T128);
627
  }
628
 
629
  ed2 = 33;
630
  if (__unsigned_compare_gt_128 (CY, CA))
631
    ed2++;
632
 
633
  T128 = power10_table_128[ed2];
634
  __mul_128x128_to_256 (CA4, CA, T128);
635
 
636
  ed2 += estimate_decimal_digits[bin_index];
637
  CQ.w[0] = CQ.w[1] = 0;
638
  diff_expon = diff_expon - ed2;
639
 
640
} else {
641
  // get CQ = CX/CY
642
  __div_128_by_128 (&CQ, &CR, CX, CY);
643
 
644
  if (!CR.w[1] && !CR.w[0]) {
645
    get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
646
                pfpsf);
647
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
648
    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
649
#endif
650
    BID_RETURN (res);
651
  }
652
  // get number of decimal digits in CQ
653
  // 2^64
654
  f64.i = 0x5f800000;
655
  fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
656
  // binary expon. of CQ
657
  bin_expon = (fx.i - 0x3f800000) >> 23;
658
 
659
  digits_q = estimate_decimal_digits[bin_expon];
660
  TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
661
  TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
662
  if (__unsigned_compare_ge_128 (CQ, TP128))
663
    digits_q++;
664
 
665
  ed2 = 34 - digits_q;
666
  T128.w[0] = power10_table_128[ed2].w[0];
667
  T128.w[1] = power10_table_128[ed2].w[1];
668
  __mul_128x128_to_256 (CA4, CR, T128);
669
  diff_expon = diff_expon - ed2;
670
  __mul_128x128_low (CQ, CQ, T128);
671
 
672
}
673
 
674
__div_256_by_128 (&CQ, &CA4, CY);
675
 
676
 
677
#ifdef SET_STATUS_FLAGS
678
  if (CA4.w[0] || CA4.w[1]) {
679
    // set status flags
680
    __set_status_flags (pfpsf, INEXACT_EXCEPTION);
681
  }
682
#ifndef LEAVE_TRAILING_ZEROS
683
  else
684
#endif
685
#else
686
#ifndef LEAVE_TRAILING_ZEROS
687
  if (!CA4.w[0] && !CA4.w[1])
688
#endif
689
#endif
690
#ifndef LEAVE_TRAILING_ZEROS
691
    // check whether result is exact
692
  {
693
    // check whether CX, CY are short
694
    if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
695
      i = (int) CY.w[0] - 1;
696
      j = (int) CX.w[0] - 1;
697
      // difference in powers of 2 factors for Y and X
698
      nzeros = ed2 - factors[i][0] + factors[j][0];
699
      // difference in powers of 5 factors
700
      d5 = ed2 - factors[i][1] + factors[j][1];
701
      if (d5 < nzeros)
702
        nzeros = d5;
703
      // get P*(2^M[extra_digits])/10^extra_digits
704
      __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
705
      //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
706
 
707
      // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
708
      amount = recip_scale[nzeros];
709
      __shr_128_long (CQ, Qh, amount);
710
 
711
      diff_expon += nzeros;
712
    } else {
713
      // decompose Q as Qh*10^17 + Ql
714
      //T128 = reciprocals10_128[17];
715
      T128.w[0] = 0x44909befeb9fad49ull;
716
      T128.w[1] = 0x000b877aa3236a4bull;
717
      __mul_128x128_to_256 (P256, CQ, T128);
718
      //amount = recip_scale[17];
719
      Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
720
      Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
721
 
722
      if (!Q_low) {
723
        diff_expon += 17;
724
 
725
        tdigit[0] = Q_high & 0x3ffffff;
726
        tdigit[1] = 0;
727
        QX = Q_high >> 26;
728
        QX32 = QX;
729
        nzeros = 0;
730
 
731
        for (j = 0; QX32; j++, QX32 >>= 7) {
732
          k = (QX32 & 127);
733
          tdigit[0] += convert_table[j][k][0];
734
          tdigit[1] += convert_table[j][k][1];
735
          if (tdigit[0] >= 100000000) {
736
            tdigit[0] -= 100000000;
737
            tdigit[1]++;
738
          }
739
        }
740
 
741
 
742
        if (tdigit[1] >= 100000000) {
743
          tdigit[1] -= 100000000;
744
          if (tdigit[1] >= 100000000)
745
            tdigit[1] -= 100000000;
746
        }
747
 
748
        digit = tdigit[0];
749
        if (!digit && !tdigit[1])
750
          nzeros += 16;
751
        else {
752
          if (!digit) {
753
            nzeros += 8;
754
            digit = tdigit[1];
755
          }
756
          // decompose digit
757
          PD = (UINT64) digit *0x068DB8BBull;
758
          digit_h = (UINT32) (PD >> 40);
759
          digit_low = digit - digit_h * 10000;
760
 
761
          if (!digit_low)
762
            nzeros += 4;
763
          else
764
            digit_h = digit_low;
765
 
766
          if (!(digit_h & 1))
767
            nzeros +=
768
              3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
769
                            (digit_h & 7));
770
        }
771
 
772
        if (nzeros) {
773
          __mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
774
 
775
          // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
776
          amount = short_recip_scale[nzeros];
777
          CQ.w[0] = CQ.w[1] >> amount;
778
        } else
779
          CQ.w[0] = Q_high;
780
        CQ.w[1] = 0;
781
 
782
        diff_expon += nzeros;
783
      } else {
784
        tdigit[0] = Q_low & 0x3ffffff;
785
        tdigit[1] = 0;
786
        QX = Q_low >> 26;
787
        QX32 = QX;
788
        nzeros = 0;
789
 
790
        for (j = 0; QX32; j++, QX32 >>= 7) {
791
          k = (QX32 & 127);
792
          tdigit[0] += convert_table[j][k][0];
793
          tdigit[1] += convert_table[j][k][1];
794
          if (tdigit[0] >= 100000000) {
795
            tdigit[0] -= 100000000;
796
            tdigit[1]++;
797
          }
798
        }
799
 
800
        if (tdigit[1] >= 100000000) {
801
          tdigit[1] -= 100000000;
802
          if (tdigit[1] >= 100000000)
803
            tdigit[1] -= 100000000;
804
        }
805
 
806
        digit = tdigit[0];
807
        if (!digit && !tdigit[1])
808
          nzeros += 16;
809
        else {
810
          if (!digit) {
811
            nzeros += 8;
812
            digit = tdigit[1];
813
          }
814
          // decompose digit
815
          PD = (UINT64) digit *0x068DB8BBull;
816
          digit_h = (UINT32) (PD >> 40);
817
          digit_low = digit - digit_h * 10000;
818
 
819
          if (!digit_low)
820
            nzeros += 4;
821
          else
822
            digit_h = digit_low;
823
 
824
          if (!(digit_h & 1))
825
            nzeros +=
826
              3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
827
                            (digit_h & 7));
828
        }
829
 
830
        if (nzeros) {
831
          // get P*(2^M[extra_digits])/10^extra_digits
832
          __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
833
 
834
          // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
835
          amount = recip_scale[nzeros];
836
          __shr_128 (CQ, Qh, amount);
837
        }
838
        diff_expon += nzeros;
839
 
840
      }
841
    }
842
    get_BID128(&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,pfpsf);
843
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
844
    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
845
#endif
846
    BID_RETURN (res);
847
  }
848
#endif
849
 
850
if (diff_expon >= 0) {
851
#ifdef IEEE_ROUND_NEAREST
852
  // rounding
853
  // 2*CA4 - CY
854
  CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
855
  CA4r.w[0] = CA4.w[0] + CA4.w[0];
856
  __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
857
  CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
858
 
859
  D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
860
  carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
861
 
862
  CQ.w[0] += carry64;
863
  if (CQ.w[0] < carry64)
864
    CQ.w[1]++;
865
#else
866
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
867
  // rounding
868
  // 2*CA4 - CY
869
  CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
870
  CA4r.w[0] = CA4.w[0] + CA4.w[0];
871
  __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
872
  CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
873
 
874
  D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
875
  carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
876
 
877
  CQ.w[0] += carry64;
878
  if (CQ.w[0] < carry64)
879
    CQ.w[1]++;
880
#else
881
  rmode = rnd_mode;
882
  if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
883
    rmode = 3 - rmode;
884
  switch (rmode) {
885
  case ROUNDING_TO_NEAREST:     // round to nearest code
886
    // rounding
887
    // 2*CA4 - CY
888
    CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
889
    CA4r.w[0] = CA4.w[0] + CA4.w[0];
890
    __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
891
    CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
892
    D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
893
    carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
894
    CQ.w[0] += carry64;
895
    if (CQ.w[0] < carry64)
896
      CQ.w[1]++;
897
    break;
898
  case ROUNDING_TIES_AWAY:
899
    // rounding
900
    // 2*CA4 - CY
901
    CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
902
    CA4r.w[0] = CA4.w[0] + CA4.w[0];
903
    __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
904
    CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
905
    D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
906
    carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
907
    CQ.w[0] += carry64;
908
    if (CQ.w[0] < carry64)
909
      CQ.w[1]++;
910
    break;
911
  case ROUNDING_DOWN:
912
  case ROUNDING_TO_ZERO:
913
    break;
914
  default:      // rounding up
915
    CQ.w[0]++;
916
    if (!CQ.w[0])
917
      CQ.w[1]++;
918
    break;
919
  }
920
#endif
921
#endif
922
 
923
} else {
924
#ifdef SET_STATUS_FLAGS
925
  if (CA4.w[0] || CA4.w[1]) {
926
    // set status flags
927
    __set_status_flags (pfpsf, INEXACT_EXCEPTION);
928
  }
929
#endif
930
  handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
931
                     CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
932
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
933
  (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
934
#endif
935
  BID_RETURN (res);
936
 
937
}
938
 
939
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
940
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
941
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
942
#endif
943
BID_RETURN (res);
944
}
945
 
946
 
947
BID128_FUNCTION_ARGTYPE1_ARG128 (bid128dq_div, UINT64, x, y)
948
     UINT256 CA4, CA4r, P256;
949
     UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
950
     UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, valid_y,
951
       PD;
952
     int_float fx, fy, f64;
953
     UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
954
     int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
955
       digits_q, amount;
956
     int nzeros, i, j, k, d5;
957
     unsigned rmode;
958
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
959
     fexcept_t binaryflags = 0;
960
#endif
961
 
962
valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
963
 
964
        // unpack arguments, check for NaN or Infinity
965
CX.w[1] = 0;
966
if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], x)) {
967
#ifdef SET_STATUS_FLAGS
968
if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64)      // y is sNaN
969
  __set_status_flags (pfpsf, INVALID_EXCEPTION);
970
#endif
971
 
972
    // test if x is NaN
973
if ((x & NAN_MASK64) == NAN_MASK64) {
974
#ifdef SET_STATUS_FLAGS
975
  if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
976
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
977
#endif
978
  res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
979
  __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
980
  res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
981
  BID_RETURN (res);
982
}
983
           // x is Infinity?
984
if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
985
  // check if y is Inf.
986
  if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
987
    // return NaN 
988
  {
989
#ifdef SET_STATUS_FLAGS
990
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
991
#endif
992
    res.w[1] = 0x7c00000000000000ull;
993
    res.w[0] = 0;
994
    BID_RETURN (res);
995
  }
996
  if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
997
  // otherwise return +/-Inf
998
  res.w[1] =
999
    ((x ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
1000
  res.w[0] = 0;
1001
  BID_RETURN (res);
1002
  }
1003
}
1004
           // x is 0
1005
if ((y.w[1] & INFINITY_MASK64) != INFINITY_MASK64) {
1006
  if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) {
1007
#ifdef SET_STATUS_FLAGS
1008
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
1009
#endif
1010
    // x=y=0, return NaN
1011
    res.w[1] = 0x7c00000000000000ull;
1012
    res.w[0] = 0;
1013
    BID_RETURN (res);
1014
  }
1015
  // return 0
1016
  res.w[1] = (x ^ y.w[1]) & 0x8000000000000000ull;
1017
  exponent_x = exponent_x - exponent_y + (DECIMAL_EXPONENT_BIAS_128<<1) - DECIMAL_EXPONENT_BIAS;
1018
  if (exponent_x > DECIMAL_MAX_EXPON_128)
1019
    exponent_x = DECIMAL_MAX_EXPON_128;
1020
  else if (exponent_x < 0)
1021
    exponent_x = 0;
1022
  res.w[1] |= (((UINT64) exponent_x) << 49);
1023
  res.w[0] = 0;
1024
  BID_RETURN (res);
1025
}
1026
}
1027
exponent_x += (DECIMAL_EXPONENT_BIAS_128 - DECIMAL_EXPONENT_BIAS);
1028
 
1029
if (!valid_y) {
1030
  // y is Inf. or NaN
1031
 
1032
  // test if y is NaN
1033
  if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
1034
#ifdef SET_STATUS_FLAGS
1035
    if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull)      // sNaN
1036
      __set_status_flags (pfpsf, INVALID_EXCEPTION);
1037
#endif
1038
    res.w[1] = CY.w[1] & QUIET_MASK64;
1039
    res.w[0] = CY.w[0];
1040
    BID_RETURN (res);
1041
  }
1042
  // y is Infinity?
1043
  if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
1044
    // return +/-0
1045
    res.w[1] = sign_x ^ sign_y;
1046
    res.w[0] = 0;
1047
    BID_RETURN (res);
1048
  }
1049
  // y is 0, return +/-Inf
1050
  res.w[1] =
1051
    ((x ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
1052
  res.w[0] = 0;
1053
#ifdef SET_STATUS_FLAGS
1054
  __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
1055
#endif
1056
  BID_RETURN (res);
1057
}
1058
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1059
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
1060
#endif
1061
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
1062
 
1063
if (__unsigned_compare_gt_128 (CY, CX)) {
1064
  // CX < CY
1065
 
1066
  // 2^64
1067
  f64.i = 0x5f800000;
1068
 
1069
  // fx ~ CX,   fy ~ CY
1070
  fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
1071
  fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
1072
  // expon_cy - expon_cx
1073
  bin_index = (fy.i - fx.i) >> 23;
1074
 
1075
  if (CX.w[1]) {
1076
    T = power10_index_binexp_128[bin_index].w[0];
1077
    __mul_64x128_short (CA, T, CX);
1078
  } else {
1079
    T128 = power10_index_binexp_128[bin_index];
1080
    __mul_64x128_short (CA, CX.w[0], T128);
1081
  }
1082
 
1083
  ed2 = 33;
1084
  if (__unsigned_compare_gt_128 (CY, CA))
1085
    ed2++;
1086
 
1087
  T128 = power10_table_128[ed2];
1088
  __mul_128x128_to_256 (CA4, CA, T128);
1089
 
1090
  ed2 += estimate_decimal_digits[bin_index];
1091
  CQ.w[0] = CQ.w[1] = 0;
1092
  diff_expon = diff_expon - ed2;
1093
 
1094
} else {
1095
  // get CQ = CX/CY
1096
  __div_128_by_128 (&CQ, &CR, CX, CY);
1097
 
1098
  if (!CR.w[1] && !CR.w[0]) {
1099
    get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
1100
                pfpsf);
1101
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1102
    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1103
#endif
1104
    BID_RETURN (res);
1105
  }
1106
  // get number of decimal digits in CQ
1107
  // 2^64
1108
  f64.i = 0x5f800000;
1109
  fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
1110
  // binary expon. of CQ
1111
  bin_expon = (fx.i - 0x3f800000) >> 23;
1112
 
1113
  digits_q = estimate_decimal_digits[bin_expon];
1114
  TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
1115
  TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
1116
  if (__unsigned_compare_ge_128 (CQ, TP128))
1117
    digits_q++;
1118
 
1119
  ed2 = 34 - digits_q;
1120
  T128.w[0] = power10_table_128[ed2].w[0];
1121
  T128.w[1] = power10_table_128[ed2].w[1];
1122
  __mul_128x128_to_256 (CA4, CR, T128);
1123
  diff_expon = diff_expon - ed2;
1124
  __mul_128x128_low (CQ, CQ, T128);
1125
 
1126
}
1127
 
1128
__div_256_by_128 (&CQ, &CA4, CY);
1129
 
1130
#ifdef SET_STATUS_FLAGS
1131
  if (CA4.w[0] || CA4.w[1]) {
1132
    // set status flags
1133
    __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1134
  }
1135
#ifndef LEAVE_TRAILING_ZEROS
1136
  else
1137
#endif
1138
#else
1139
#ifndef LEAVE_TRAILING_ZEROS
1140
  if (!CA4.w[0] && !CA4.w[1])
1141
#endif
1142
#endif
1143
#ifndef LEAVE_TRAILING_ZEROS
1144
    // check whether result is exact
1145
  {
1146
    //printf("ed2=%d,nz=%d,a=%d,CQ="LX16","LX16", RH="LX16", RL="LX16"\n",ed2,nzeros,amount,CQ.w[1],CQ.w[0],reciprocals10_128[nzeros].w[1],reciprocals10_128[nzeros].w[0]);fflush(stdout);
1147
    // check whether CX, CY are short
1148
    if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
1149
      i = (int) CY.w[0] - 1;
1150
      j = (int) CX.w[0] - 1;
1151
      // difference in powers of 2 factors for Y and X
1152
      nzeros = ed2 - factors[i][0] + factors[j][0];
1153
      // difference in powers of 5 factors
1154
      d5 = ed2 - factors[i][1] + factors[j][1];
1155
      if (d5 < nzeros)
1156
        nzeros = d5;
1157
      // get P*(2^M[extra_digits])/10^extra_digits
1158
      __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
1159
      //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
1160
 
1161
      // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1162
      amount = recip_scale[nzeros];
1163
      __shr_128_long (CQ, Qh, amount);
1164
 
1165
      diff_expon += nzeros;
1166
    } else {
1167
      // decompose Q as Qh*10^17 + Ql
1168
      //T128 = reciprocals10_128[17];
1169
      T128.w[0] = 0x44909befeb9fad49ull;
1170
      T128.w[1] = 0x000b877aa3236a4bull;
1171
      __mul_128x128_to_256 (P256, CQ, T128);
1172
      //amount = recip_scale[17];
1173
      Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
1174
      Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
1175
 
1176
      if (!Q_low) {
1177
        diff_expon += 17;
1178
 
1179
        tdigit[0] = Q_high & 0x3ffffff;
1180
        tdigit[1] = 0;
1181
        QX = Q_high >> 26;
1182
        QX32 = QX;
1183
        nzeros = 0;
1184
 
1185
        for (j = 0; QX32; j++, QX32 >>= 7) {
1186
          k = (QX32 & 127);
1187
          tdigit[0] += convert_table[j][k][0];
1188
          tdigit[1] += convert_table[j][k][1];
1189
          if (tdigit[0] >= 100000000) {
1190
            tdigit[0] -= 100000000;
1191
            tdigit[1]++;
1192
          }
1193
        }
1194
 
1195
 
1196
        if (tdigit[1] >= 100000000) {
1197
          tdigit[1] -= 100000000;
1198
          if (tdigit[1] >= 100000000)
1199
            tdigit[1] -= 100000000;
1200
        }
1201
 
1202
        digit = tdigit[0];
1203
        if (!digit && !tdigit[1])
1204
          nzeros += 16;
1205
        else {
1206
          if (!digit) {
1207
            nzeros += 8;
1208
            digit = tdigit[1];
1209
          }
1210
          // decompose digit
1211
          PD = (UINT64) digit *0x068DB8BBull;
1212
          digit_h = (UINT32) (PD >> 40);
1213
          //printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
1214
          digit_low = digit - digit_h * 10000;
1215
 
1216
          if (!digit_low)
1217
            nzeros += 4;
1218
          else
1219
            digit_h = digit_low;
1220
 
1221
          if (!(digit_h & 1))
1222
            nzeros +=
1223
              3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
1224
                            (digit_h & 7));
1225
        }
1226
 
1227
        if (nzeros) {
1228
          __mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
1229
 
1230
          // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
1231
          amount = short_recip_scale[nzeros];
1232
          CQ.w[0] = CQ.w[1] >> amount;
1233
        } else
1234
          CQ.w[0] = Q_high;
1235
        CQ.w[1] = 0;
1236
 
1237
        diff_expon += nzeros;
1238
      } else {
1239
        tdigit[0] = Q_low & 0x3ffffff;
1240
        tdigit[1] = 0;
1241
        QX = Q_low >> 26;
1242
        QX32 = QX;
1243
        nzeros = 0;
1244
 
1245
        for (j = 0; QX32; j++, QX32 >>= 7) {
1246
          k = (QX32 & 127);
1247
          tdigit[0] += convert_table[j][k][0];
1248
          tdigit[1] += convert_table[j][k][1];
1249
          if (tdigit[0] >= 100000000) {
1250
            tdigit[0] -= 100000000;
1251
            tdigit[1]++;
1252
          }
1253
        }
1254
 
1255
        if (tdigit[1] >= 100000000) {
1256
          tdigit[1] -= 100000000;
1257
          if (tdigit[1] >= 100000000)
1258
            tdigit[1] -= 100000000;
1259
        }
1260
 
1261
        digit = tdigit[0];
1262
        if (!digit && !tdigit[1])
1263
          nzeros += 16;
1264
        else {
1265
          if (!digit) {
1266
            nzeros += 8;
1267
            digit = tdigit[1];
1268
          }
1269
          // decompose digit
1270
          PD = (UINT64) digit *0x068DB8BBull;
1271
          digit_h = (UINT32) (PD >> 40);
1272
          //printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
1273
          digit_low = digit - digit_h * 10000;
1274
 
1275
          if (!digit_low)
1276
            nzeros += 4;
1277
          else
1278
            digit_h = digit_low;
1279
 
1280
          if (!(digit_h & 1))
1281
            nzeros +=
1282
              3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
1283
                            (digit_h & 7));
1284
        }
1285
 
1286
        if (nzeros) {
1287
          // get P*(2^M[extra_digits])/10^extra_digits
1288
          __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
1289
 
1290
          // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1291
          amount = recip_scale[nzeros];
1292
          __shr_128 (CQ, Qh, amount);
1293
        }
1294
        diff_expon += nzeros;
1295
 
1296
      }
1297
    }
1298
    get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
1299
                pfpsf);
1300
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1301
    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1302
#endif
1303
    BID_RETURN (res);
1304
  }
1305
#endif
1306
 
1307
if (diff_expon >= 0) {
1308
#ifdef IEEE_ROUND_NEAREST
1309
  // rounding
1310
  // 2*CA4 - CY
1311
  CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1312
  CA4r.w[0] = CA4.w[0] + CA4.w[0];
1313
  __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1314
  CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1315
 
1316
  D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1317
  carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1318
 
1319
  CQ.w[0] += carry64;
1320
  if (CQ.w[0] < carry64)
1321
    CQ.w[1]++;
1322
#else
1323
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
1324
  // rounding
1325
  // 2*CA4 - CY
1326
  CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1327
  CA4r.w[0] = CA4.w[0] + CA4.w[0];
1328
  __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1329
  CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1330
 
1331
  D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1332
  carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1333
 
1334
  CQ.w[0] += carry64;
1335
  if (CQ.w[0] < carry64)
1336
    CQ.w[1]++;
1337
#else
1338
  rmode = rnd_mode;
1339
  if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
1340
    rmode = 3 - rmode;
1341
  switch (rmode) {
1342
  case ROUNDING_TO_NEAREST:     // round to nearest code
1343
    // rounding
1344
    // 2*CA4 - CY
1345
    CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1346
    CA4r.w[0] = CA4.w[0] + CA4.w[0];
1347
    __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1348
    CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1349
    D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1350
    carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1351
    CQ.w[0] += carry64;
1352
    if (CQ.w[0] < carry64)
1353
      CQ.w[1]++;
1354
    break;
1355
  case ROUNDING_TIES_AWAY:
1356
    // rounding
1357
    // 2*CA4 - CY
1358
    CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1359
    CA4r.w[0] = CA4.w[0] + CA4.w[0];
1360
    __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1361
    CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1362
    D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1363
    carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1364
    CQ.w[0] += carry64;
1365
    if (CQ.w[0] < carry64)
1366
      CQ.w[1]++;
1367
    break;
1368
  case ROUNDING_DOWN:
1369
  case ROUNDING_TO_ZERO:
1370
    break;
1371
  default:      // rounding up
1372
    CQ.w[0]++;
1373
    if (!CQ.w[0])
1374
      CQ.w[1]++;
1375
    break;
1376
  }
1377
#endif
1378
#endif
1379
 
1380
} else {
1381
#ifdef SET_STATUS_FLAGS
1382
  if (CA4.w[0] || CA4.w[1]) {
1383
    // set status flags
1384
    __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1385
  }
1386
#endif
1387
  handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
1388
                     CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
1389
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1390
  (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1391
#endif
1392
  BID_RETURN (res);
1393
}
1394
 
1395
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
1396
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1397
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1398
#endif
1399
BID_RETURN (res);
1400
 
1401
}
1402
 
1403
 
1404
BID128_FUNCTION_ARG128_ARGTYPE2 (bid128qd_div, x, UINT64, y)
1405
     UINT256 CA4, CA4r, P256;
1406
     UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
1407
     UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, PD,
1408
       valid_y;
1409
     int_float fx, fy, f64;
1410
     UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
1411
     int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
1412
       digits_q, amount;
1413
     int nzeros, i, j, k, d5, rmode;
1414
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1415
     fexcept_t binaryflags = 0;
1416
#endif
1417
 
1418
 
1419
valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], y);
1420
        // unpack arguments, check for NaN or Infinity
1421
if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
1422
    // test if x is NaN
1423
if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
1424
#ifdef SET_STATUS_FLAGS
1425
  if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull ||      // sNaN
1426
      (y & 0x7e00000000000000ull) == 0x7e00000000000000ull)
1427
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
1428
#endif
1429
  res.w[1] = (CX.w[1]) & QUIET_MASK64;
1430
  res.w[0] = CX.w[0];
1431
  BID_RETURN (res);
1432
}
1433
    // x is Infinity?
1434
if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
1435
  // check if y is Inf. 
1436
  if (((y & 0x7c00000000000000ull) == 0x7800000000000000ull))
1437
    // return NaN 
1438
  {
1439
#ifdef SET_STATUS_FLAGS
1440
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
1441
#endif
1442
    res.w[1] = 0x7c00000000000000ull;
1443
    res.w[0] = 0;
1444
    BID_RETURN (res);
1445
  }
1446
  // y is NaN?
1447
  if (((y & 0x7c00000000000000ull) != 0x7c00000000000000ull))
1448
    // return NaN 
1449
  {
1450
    // return +/-Inf
1451
    res.w[1] = ((x.w[1] ^ y) & 0x8000000000000000ull) |
1452
      0x7800000000000000ull;
1453
    res.w[0] = 0;
1454
    BID_RETURN (res);
1455
  }
1456
}
1457
    // x is 0
1458
if ((y & 0x7800000000000000ull) < 0x7800000000000000ull) {
1459
        if (!CY.w[0]) {
1460
#ifdef SET_STATUS_FLAGS
1461
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
1462
#endif
1463
    // x=y=0, return NaN
1464
    res.w[1] = 0x7c00000000000000ull;
1465
    res.w[0] = 0;
1466
    BID_RETURN (res);
1467
  }
1468
  // return 0
1469
  res.w[1] = (x.w[1] ^ y) & 0x8000000000000000ull;
1470
  exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
1471
  if (exponent_x > DECIMAL_MAX_EXPON_128)
1472
    exponent_x = DECIMAL_MAX_EXPON_128;
1473
  else if (exponent_x < 0)
1474
    exponent_x = 0;
1475
  res.w[1] |= (((UINT64) exponent_x) << 49);
1476
  res.w[0] = 0;
1477
  BID_RETURN (res);
1478
}
1479
}
1480
CY.w[1] = 0;
1481
if (!valid_y) {
1482
  // y is Inf. or NaN
1483
 
1484
  // test if y is NaN
1485
  if ((y & NAN_MASK64) == NAN_MASK64) {
1486
#ifdef SET_STATUS_FLAGS
1487
    if ((y & SNAN_MASK64) == SNAN_MASK64)       // sNaN
1488
      __set_status_flags (pfpsf, INVALID_EXCEPTION);
1489
#endif
1490
  res.w[0] = (CY.w[0] & 0x0003ffffffffffffull);
1491
  __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
1492
  res.w[1] |= ((CY.w[0]) & 0xfc00000000000000ull);
1493
    BID_RETURN (res);
1494
  }
1495
  // y is Infinity?
1496
  if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
1497
    // return +/-0
1498
    res.w[1] = ((x.w[1] ^ y) & 0x8000000000000000ull);
1499
    res.w[0] = 0;
1500
    BID_RETURN (res);
1501
  }
1502
  // y is 0
1503
#ifdef SET_STATUS_FLAGS
1504
  __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
1505
#endif
1506
  res.w[1] = (sign_x ^ sign_y) | INFINITY_MASK64;
1507
  res.w[0] = 0;
1508
  BID_RETURN (res);
1509
}
1510
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1511
(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
1512
#endif
1513
diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
1514
 
1515
if (__unsigned_compare_gt_128 (CY, CX)) {
1516
  // CX < CY
1517
 
1518
  // 2^64
1519
  f64.i = 0x5f800000;
1520
 
1521
  // fx ~ CX,   fy ~ CY
1522
  fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
1523
  fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
1524
  // expon_cy - expon_cx
1525
  bin_index = (fy.i - fx.i) >> 23;
1526
 
1527
  if (CX.w[1]) {
1528
    T = power10_index_binexp_128[bin_index].w[0];
1529
    __mul_64x128_short (CA, T, CX);
1530
  } else {
1531
    T128 = power10_index_binexp_128[bin_index];
1532
    __mul_64x128_short (CA, CX.w[0], T128);
1533
  }
1534
 
1535
  ed2 = 33;
1536
  if (__unsigned_compare_gt_128 (CY, CA))
1537
    ed2++;
1538
 
1539
  T128 = power10_table_128[ed2];
1540
  __mul_128x128_to_256 (CA4, CA, T128);
1541
 
1542
  ed2 += estimate_decimal_digits[bin_index];
1543
  CQ.w[0] = CQ.w[1] = 0;
1544
  diff_expon = diff_expon - ed2;
1545
 
1546
} else {
1547
  // get CQ = CX/CY
1548
  __div_128_by_128 (&CQ, &CR, CX, CY);
1549
 
1550
  if (!CR.w[1] && !CR.w[0]) {
1551
    get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
1552
                pfpsf);
1553
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1554
    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1555
#endif
1556
    BID_RETURN (res);
1557
  }
1558
  // get number of decimal digits in CQ
1559
  // 2^64
1560
  f64.i = 0x5f800000;
1561
  fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
1562
  // binary expon. of CQ
1563
  bin_expon = (fx.i - 0x3f800000) >> 23;
1564
 
1565
  digits_q = estimate_decimal_digits[bin_expon];
1566
  TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
1567
  TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
1568
  if (__unsigned_compare_ge_128 (CQ, TP128))
1569
    digits_q++;
1570
 
1571
  ed2 = 34 - digits_q;
1572
  T128.w[0] = power10_table_128[ed2].w[0];
1573
  T128.w[1] = power10_table_128[ed2].w[1];
1574
  __mul_128x128_to_256 (CA4, CR, T128);
1575
  diff_expon = diff_expon - ed2;
1576
  __mul_128x128_low (CQ, CQ, T128);
1577
 
1578
}
1579
 
1580
__div_256_by_128 (&CQ, &CA4, CY);
1581
 
1582
 
1583
#ifdef SET_STATUS_FLAGS
1584
  if (CA4.w[0] || CA4.w[1]) {
1585
    // set status flags
1586
    __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1587
  }
1588
#ifndef LEAVE_TRAILING_ZEROS
1589
  else
1590
#endif
1591
#else
1592
#ifndef LEAVE_TRAILING_ZEROS
1593
  if (!CA4.w[0] && !CA4.w[1])
1594
#endif
1595
#endif
1596
#ifndef LEAVE_TRAILING_ZEROS
1597
    // check whether result is exact
1598
  {
1599
    // check whether CX, CY are short
1600
    if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
1601
      i = (int) CY.w[0] - 1;
1602
      j = (int) CX.w[0] - 1;
1603
      // difference in powers of 2 factors for Y and X
1604
      nzeros = ed2 - factors[i][0] + factors[j][0];
1605
      // difference in powers of 5 factors
1606
      d5 = ed2 - factors[i][1] + factors[j][1];
1607
      if (d5 < nzeros)
1608
        nzeros = d5;
1609
      // get P*(2^M[extra_digits])/10^extra_digits
1610
      __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
1611
      //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
1612
 
1613
      // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1614
      amount = recip_scale[nzeros];
1615
      __shr_128_long (CQ, Qh, amount);
1616
 
1617
      diff_expon += nzeros;
1618
    } else {
1619
      // decompose Q as Qh*10^17 + Ql
1620
      //T128 = reciprocals10_128[17];
1621
      T128.w[0] = 0x44909befeb9fad49ull;
1622
      T128.w[1] = 0x000b877aa3236a4bull;
1623
      __mul_128x128_to_256 (P256, CQ, T128);
1624
      //amount = recip_scale[17];
1625
      Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
1626
      Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
1627
 
1628
      if (!Q_low) {
1629
        diff_expon += 17;
1630
 
1631
        tdigit[0] = Q_high & 0x3ffffff;
1632
        tdigit[1] = 0;
1633
        QX = Q_high >> 26;
1634
        QX32 = QX;
1635
        nzeros = 0;
1636
 
1637
        for (j = 0; QX32; j++, QX32 >>= 7) {
1638
          k = (QX32 & 127);
1639
          tdigit[0] += convert_table[j][k][0];
1640
          tdigit[1] += convert_table[j][k][1];
1641
          if (tdigit[0] >= 100000000) {
1642
            tdigit[0] -= 100000000;
1643
            tdigit[1]++;
1644
          }
1645
        }
1646
 
1647
 
1648
        if (tdigit[1] >= 100000000) {
1649
          tdigit[1] -= 100000000;
1650
          if (tdigit[1] >= 100000000)
1651
            tdigit[1] -= 100000000;
1652
        }
1653
 
1654
        digit = tdigit[0];
1655
        if (!digit && !tdigit[1])
1656
          nzeros += 16;
1657
        else {
1658
          if (!digit) {
1659
            nzeros += 8;
1660
            digit = tdigit[1];
1661
          }
1662
          // decompose digit
1663
          PD = (UINT64) digit *0x068DB8BBull;
1664
          digit_h = (UINT32) (PD >> 40);
1665
          digit_low = digit - digit_h * 10000;
1666
 
1667
          if (!digit_low)
1668
            nzeros += 4;
1669
          else
1670
            digit_h = digit_low;
1671
 
1672
          if (!(digit_h & 1))
1673
            nzeros +=
1674
              3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
1675
                            (digit_h & 7));
1676
        }
1677
 
1678
        if (nzeros) {
1679
          __mul_64x64_to_128 (CQ, Q_high, reciprocals10_64[nzeros]);
1680
 
1681
          // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
1682
          amount = short_recip_scale[nzeros];
1683
          CQ.w[0] = CQ.w[1] >> amount;
1684
        } else
1685
          CQ.w[0] = Q_high;
1686
        CQ.w[1] = 0;
1687
 
1688
        diff_expon += nzeros;
1689
      } else {
1690
        tdigit[0] = Q_low & 0x3ffffff;
1691
        tdigit[1] = 0;
1692
        QX = Q_low >> 26;
1693
        QX32 = QX;
1694
        nzeros = 0;
1695
 
1696
        for (j = 0; QX32; j++, QX32 >>= 7) {
1697
          k = (QX32 & 127);
1698
          tdigit[0] += convert_table[j][k][0];
1699
          tdigit[1] += convert_table[j][k][1];
1700
          if (tdigit[0] >= 100000000) {
1701
            tdigit[0] -= 100000000;
1702
            tdigit[1]++;
1703
          }
1704
        }
1705
 
1706
        if (tdigit[1] >= 100000000) {
1707
          tdigit[1] -= 100000000;
1708
          if (tdigit[1] >= 100000000)
1709
            tdigit[1] -= 100000000;
1710
        }
1711
 
1712
        digit = tdigit[0];
1713
        if (!digit && !tdigit[1])
1714
          nzeros += 16;
1715
        else {
1716
          if (!digit) {
1717
            nzeros += 8;
1718
            digit = tdigit[1];
1719
          }
1720
          // decompose digit
1721
          PD = (UINT64) digit *0x068DB8BBull;
1722
          digit_h = (UINT32) (PD >> 40);
1723
          digit_low = digit - digit_h * 10000;
1724
 
1725
          if (!digit_low)
1726
            nzeros += 4;
1727
          else
1728
            digit_h = digit_low;
1729
 
1730
          if (!(digit_h & 1))
1731
            nzeros +=
1732
              3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
1733
                            (digit_h & 7));
1734
        }
1735
 
1736
        if (nzeros) {
1737
          // get P*(2^M[extra_digits])/10^extra_digits
1738
          __mul_128x128_full (Qh, Ql, CQ, reciprocals10_128[nzeros]);
1739
 
1740
          // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1741
          amount = recip_scale[nzeros];
1742
          __shr_128 (CQ, Qh, amount);
1743
        }
1744
        diff_expon += nzeros;
1745
 
1746
      }
1747
    }
1748
    get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,pfpsf);
1749
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1750
    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1751
#endif
1752
    BID_RETURN (res);
1753
  }
1754
#endif
1755
 
1756
if (diff_expon >= 0) {
1757
#ifdef IEEE_ROUND_NEAREST
1758
  // rounding
1759
  // 2*CA4 - CY
1760
  CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1761
  CA4r.w[0] = CA4.w[0] + CA4.w[0];
1762
  __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1763
  CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1764
 
1765
  D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1766
  carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1767
 
1768
  CQ.w[0] += carry64;
1769
  if (CQ.w[0] < carry64)
1770
    CQ.w[1]++;
1771
#else
1772
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
1773
  // rounding
1774
  // 2*CA4 - CY
1775
  CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1776
  CA4r.w[0] = CA4.w[0] + CA4.w[0];
1777
  __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1778
  CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1779
 
1780
  D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1781
  carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1782
 
1783
  CQ.w[0] += carry64;
1784
  if (CQ.w[0] < carry64)
1785
    CQ.w[1]++;
1786
#else
1787
  rmode = rnd_mode;
1788
  if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
1789
    rmode = 3 - rmode;
1790
  switch (rmode) {
1791
  case ROUNDING_TO_NEAREST:     // round to nearest code
1792
    // rounding
1793
    // 2*CA4 - CY
1794
    CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1795
    CA4r.w[0] = CA4.w[0] + CA4.w[0];
1796
    __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1797
    CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1798
    D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1799
    carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1800
    CQ.w[0] += carry64;
1801
    if (CQ.w[0] < carry64)
1802
      CQ.w[1]++;
1803
    break;
1804
  case ROUNDING_TIES_AWAY:
1805
    // rounding
1806
    // 2*CA4 - CY
1807
    CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1808
    CA4r.w[0] = CA4.w[0] + CA4.w[0];
1809
    __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1810
    CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1811
    D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1812
    carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1813
    CQ.w[0] += carry64;
1814
    if (CQ.w[0] < carry64)
1815
      CQ.w[1]++;
1816
    break;
1817
  case ROUNDING_DOWN:
1818
  case ROUNDING_TO_ZERO:
1819
    break;
1820
  default:      // rounding up
1821
    CQ.w[0]++;
1822
    if (!CQ.w[0])
1823
      CQ.w[1]++;
1824
    break;
1825
  }
1826
#endif
1827
#endif
1828
 
1829
} else {
1830
#ifdef SET_STATUS_FLAGS
1831
  if (CA4.w[0] || CA4.w[1]) {
1832
    // set status flags
1833
    __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1834
  }
1835
#endif
1836
  handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
1837
                     CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
1838
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1839
  (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1840
#endif
1841
  BID_RETURN (res);
1842
 
1843
}
1844
 
1845
get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
1846
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
1847
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1848
#endif
1849
BID_RETURN (res);
1850
 
1851
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.