1/* Copyright (C) 2007, 2009  Free Software Foundation, Inc.
2
3This file is part of GCC.
4
5GCC is free software; you can redistribute it and/or modify it under
6the terms of the GNU General Public License as published by the Free
7Software Foundation; either version 3, or (at your option) any later
8version.
9
10GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11WARRANTY; without even the implied warranty of MERCHANTABILITY or
12FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
13for more details.
14
15Under Section 7 of GPL version 3, you are granted additional
16permissions described in the GCC Runtime Library Exception, version
173.1, as published by the Free Software Foundation.
18
19You should have received a copy of the GNU General Public License and
20a copy of the GCC Runtime Library Exception along with this program;
21see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
22<http://www.gnu.org/licenses/>.  */
23
24#define  DECNUMDIGITS 34	// work with up to 34 digits
25
26#include "bid_internal.h"
27#include "bid_b2d.h"
28
29UINT32
30bid_to_bid32 (UINT32 ba) {
31  UINT32 res;
32  UINT32 sign, comb, exp;
33  UINT32 trailing;
34  UINT32 bcoeff;
35
36  sign = (ba & 0x80000000ul);
37  comb = (ba & 0x7ff00000ul) >> 20;
38  trailing = (ba & 0x000ffffful);
39
40  if ((comb & 0x780) == 0x780) {
41    ba &= 0xfff0000ul;
42    return ba;
43  } else {
44    if ((comb & 0x600) == 0x600) {	// G0..G1 = 11 -> exp is G2..G11
45      exp = (comb >> 1) & 0xff;
46      bcoeff = ((8 + (comb & 1)) << 20) | trailing;
47    } else {
48      exp = (comb >> 3) & 0xff;
49      bcoeff = ((comb & 7) << 20) | trailing;
50    }
51    if (bcoeff >= 10000000)
52      bcoeff = 0;
53    res = very_fast_get_BID32 (sign, exp, bcoeff);
54  }
55  return res;
56}
57
58UINT64
59bid_to_bid64 (UINT64 ba) {
60  UINT64 res;
61  UINT64 sign, comb, exp;
62  UINT64 trailing;
63  UINT64 bcoeff;
64
65  sign = (ba & 0x8000000000000000ull);
66  comb = (ba & 0x7ffc000000000000ull) >> 50;
67  trailing = (ba & 0x0003ffffffffffffull);
68
69  if ((comb & 0x1e00) == 0x1e00) {
70    ba &= 0xfff000000000000ULL;
71    return ba;
72  } else {
73    if ((comb & 0x1800) == 0x1800) {	// G0..G1 = 11 -> exp is G2..G11
74      exp = (comb >> 1) & 0x3ff;
75      bcoeff = ((8 + (comb & 1)) << 50) | trailing;
76    } else {
77      exp = (comb >> 3) & 0x3ff;
78      bcoeff = ((comb & 7) << 50) | trailing;
79    }
80    if (bcoeff >= 10000000000000000ull)
81      bcoeff = 0ull;
82    res = very_fast_get_BID64 (sign, exp, bcoeff);
83  }
84  return res;
85}
86
87#if DECIMAL_CALL_BY_REFERENCE
88void
89bid_to_dpd32 (UINT32 * pres, UINT32 * pba) {
90  UINT32 ba = *pba;
91#else
92UINT32
93bid_to_dpd32 (UINT32 ba) {
94#endif
95  UINT32 res;
96
97  UINT32 sign, comb, exp, trailing;
98  UINT32 b0, b1, b2;
99  UINT32 bcoeff, dcoeff;
100  UINT32 nanb = 0;
101
102  sign = (ba & 0x80000000);
103  comb = (ba & 0x7ff00000) >> 20;
104  trailing = (ba & 0xfffff);
105
106  // Detect infinity, and return canonical infinity
107  if ((comb & 0x7c0) == 0x780) {
108    res = sign | 0x78000000;
109    BID_RETURN (res);
110    // Detect NaN, and canonicalize trailing
111  } else if ((comb & 0x7c0) == 0x7c0) {
112    if (trailing > 999999)
113      trailing = 0;
114    nanb = ba & 0xfe000000;
115    exp = 0;
116    bcoeff = trailing;
117  } else {	// Normal number
118    if ((comb & 0x600) == 0x600) {	// G0..G1 = 11 -> exp is G2..G11
119      exp = (comb >> 1) & 0xff;
120      bcoeff = ((8 + (comb & 1)) << 20) | trailing;
121    } else {
122      exp = (comb >> 3) & 0xff;
123      bcoeff = ((comb & 7) << 20) | trailing;
124    }
125    // Zero the coefficient if non-canonical (>= 10^7)
126    if (bcoeff >= 10000000)
127      bcoeff = 0;
128  }
129
130  b0 = bcoeff / 1000000;
131  b1 = (bcoeff / 1000) % 1000;
132  b2 = bcoeff % 1000;
133  dcoeff = (b2d[b1] << 10) | b2d[b2];
134
135  if (b0 >= 8)	// is b0 8 or 9?
136    res =
137      sign |
138      ((0x600 | ((exp >> 6) << 7) | ((b0 & 1) << 6) | (exp & 0x3f)) <<
139       20) | dcoeff;
140  else	// else b0 is 0..7
141    res =
142      sign | ((((exp >> 6) << 9) | (b0 << 6) | (exp & 0x3f)) << 20) |
143      dcoeff;
144
145  res |= nanb;
146
147  BID_RETURN (res);
148}
149
150#if DECIMAL_CALL_BY_REFERENCE
151void
152bid_to_dpd64 (UINT64 * pres, UINT64 * pba) {
153  UINT64 ba = *pba;
154#else
155UINT64
156bid_to_dpd64 (UINT64 ba) {
157#endif
158  UINT64 res;
159
160  UINT64 sign, comb, exp;
161  UINT64 trailing;
162  UINT32 b0, b1, b2, b3, b4, b5;
163  UINT64 bcoeff;
164  UINT64 dcoeff;
165  UINT32 yhi, ylo;
166  UINT64 nanb = 0;
167
168//printf("arg bid "FMT_LLX16" \n", ba);
169  sign = (ba & 0x8000000000000000ull);
170  comb = (ba & 0x7ffc000000000000ull) >> 50;
171  trailing = (ba & 0x0003ffffffffffffull);
172
173  // Detect infinity, and return canonical infinity
174  if ((comb & 0x1f00) == 0x1e00) {
175    res = sign | 0x7800000000000000ull;
176    BID_RETURN (res);
177    // Detect NaN, and canonicalize trailing
178  } else if ((comb & 0x1e00) == 0x1e00) {
179    if (trailing > 999999999999999ull)
180      trailing = 0;
181    nanb = ba & 0xfe00000000000000ull;
182    exp = 0;
183    bcoeff = trailing;
184  } else {	// Normal number
185    if ((comb & 0x1800) == 0x1800) {	// G0..G1 = 11 -> exp is G2..G11
186      exp = (comb >> 1) & 0x3ff;
187      bcoeff = ((8 + (comb & 1)) << 50) | trailing;
188    } else {
189      exp = (comb >> 3) & 0x3ff;
190      bcoeff = ((comb & 7) << 50) | trailing;
191    }
192
193    // Zero the coefficient if it is non-canonical (>= 10^16)
194    if (bcoeff >= 10000000000000000ull)
195      bcoeff = 0;
196  }
197
198// Floor(2^61 / 10^9)
199#define D61 (2305843009ull)
200
201// Multipy the binary coefficient by ceil(2^64 / 1000), and take the upper
202// 64-bits in order to compute a division by 1000.
203
204#if 1
205  yhi =
206    ((UINT64) D61 *
207     (UINT64) (UINT32) (bcoeff >> (UINT64) 27)) >> (UINT64) 34;
208  ylo = bcoeff - 1000000000ull * yhi;
209  if (ylo >= 1000000000) {
210    ylo = ylo - 1000000000;
211    yhi = yhi + 1;
212  }
213#else
214  yhi = bcoeff / 1000000000ull;
215  ylo = bcoeff % 1000000000ull;
216#endif
217
218  // yhi = ABBBCCC ylo = DDDEEEFFF
219  b5 = ylo % 1000;	// b5 = FFF
220  b3 = ylo / 1000000;	// b3 = DDD
221  b4 = (ylo / 1000) - (1000 * b3);	// b4 = EEE
222  b2 = yhi % 1000;	// b2 = CCC
223  b0 = yhi / 1000000;	// b0 = A
224  b1 = (yhi / 1000) - (1000 * b0);	// b1 = BBB
225
226  dcoeff = b2d[b5] | b2d2[b4] | b2d3[b3] | b2d4[b2] | b2d5[b1];
227
228  if (b0 >= 8)	// is b0 8 or 9?
229    res =
230      sign |
231      ((0x1800 | ((exp >> 8) << 9) | ((b0 & 1) << 8) | (exp & 0xff)) <<
232       50) | dcoeff;
233  else	// else b0 is 0..7
234    res =
235      sign | ((((exp >> 8) << 11) | (b0 << 8) | (exp & 0xff)) << 50) |
236      dcoeff;
237
238  res |= nanb;
239
240  BID_RETURN (res);
241}
242
243#if DECIMAL_CALL_BY_REFERENCE
244void
245dpd_to_bid32 (UINT32 * pres, UINT32 * pda) {
246  UINT32 da = *pda;
247#else
248UINT32
249dpd_to_bid32 (UINT32 da) {
250#endif
251  UINT32 in = *(UINT32 *) & da;
252  UINT32 res;
253
254  UINT32 sign, comb, exp;
255  UINT32 trailing;
256  UINT32 d0 = 0, d1, d2;
257  UINT64 bcoeff;
258  UINT32 nanb = 0;
259
260  sign = (in & 0x80000000);
261  comb = (in & 0x7ff00000) >> 20;
262  trailing = (in & 0x000fffff);
263
264  if ((comb & 0x7e0) == 0x780) {	// G0..G4 = 1111 -> Inf
265    res = in & 0xf8000000;
266    BID_RETURN (res);
267  } else if ((comb & 0x7c0) == 0x7c0) {	// G0..G5 = 11111 -> NaN
268    nanb = in & 0xfe000000;
269    exp = 0;
270  } else {	// Normal number
271    if ((comb & 0x600) == 0x600) {	// G0..G1 = 11 -> d0 = 8 + G4
272      d0 = ((comb >> 6) & 1) | 8;
273      exp = ((comb & 0x180) >> 1) | (comb & 0x3f);
274    } else {
275      d0 = (comb >> 6) & 0x7;
276      exp = ((comb & 0x600) >> 3) | (comb & 0x3f);
277    }
278  }
279  d1 = d2b2[(trailing >> 10) & 0x3ff];
280  d2 = d2b[(trailing) & 0x3ff];
281
282  bcoeff = d2 + d1 + (1000000 * d0);
283  if (bcoeff < 0x800000) {
284    res = (exp << 23) | bcoeff | sign;
285  } else {
286    res = (exp << 21) | sign | 0x60000000 | (bcoeff & 0x1fffff);
287  }
288
289  res |= nanb;
290
291  BID_RETURN (res);
292}
293
294#if DECIMAL_CALL_BY_REFERENCE
295void
296dpd_to_bid64 (UINT64 * pres, UINT64 * pda) {
297  UINT64 da = *pda;
298#else
299UINT64
300dpd_to_bid64 (UINT64 da) {
301#endif
302  UINT64 in = *(UINT64 *) & da;
303  UINT64 res;
304
305  UINT64 sign, comb, exp;
306  UINT64 trailing;
307  // UINT64 d0, d1, d2, d3, d4, d5;
308
309  UINT64 d1, d2;
310  UINT32 d0, d3, d4, d5;
311  UINT64 bcoeff;
312  UINT64 nanb = 0;
313
314//printf("arg dpd "FMT_LLX16" \n", in);
315  sign = (in & 0x8000000000000000ull);
316  comb = (in & 0x7ffc000000000000ull) >> 50;
317  trailing = (in & 0x0003ffffffffffffull);
318
319  if ((comb & 0x1f00) == 0x1e00) {	// G0..G4 = 1111 -> Inf
320    res = in & 0xf800000000000000ull;
321    BID_RETURN (res);
322  } else if ((comb & 0x1f00) == 0x1f00) {	// G0..G5 = 11111 -> NaN
323    nanb = in & 0xfe00000000000000ull;
324    exp = 0;
325    d0 = 0;
326  } else {	// Normal number
327    if ((comb & 0x1800) == 0x1800) {	// G0..G1 = 11 -> d0 = 8 + G4
328      d0 = ((comb >> 8) & 1) | 8;
329      // d0 = (comb & 0x0100 ? 9 : 8);
330      exp = (comb & 0x600) >> 1;
331      // exp = (comb & 0x0400 ? 1 : 0) * 0x200 + (comb & 0x0200 ? 1 : 0) * 0x100; // exp leading bits are G2..G3
332    } else {
333      d0 = (comb >> 8) & 0x7;
334      exp = (comb & 0x1800) >> 3;
335      // exp = (comb & 0x1000 ? 1 : 0) * 0x200 + (comb & 0x0800 ? 1 : 0) * 0x100; // exp loading bits are G0..G1
336    }
337  }
338  d1 = d2b5[(trailing >> 40) & 0x3ff];
339  d2 = d2b4[(trailing >> 30) & 0x3ff];
340  d3 = d2b3[(trailing >> 20) & 0x3ff];
341  d4 = d2b2[(trailing >> 10) & 0x3ff];
342  d5 = d2b[(trailing) & 0x3ff];
343
344  bcoeff = (d5 + d4 + d3) + d2 + d1 + (1000000000000000ull * d0);
345  exp += (comb & 0xff);
346  res = very_fast_get_BID64 (sign, exp, bcoeff);
347
348  res |= nanb;
349
350  BID_RETURN (res);
351}
352
353#if DECIMAL_CALL_BY_REFERENCE
354void
355bid_to_dpd128 (UINT128 * pres, UINT128 * pba) {
356  UINT128 ba = *pba;
357#else
358UINT128
359bid_to_dpd128 (UINT128 ba) {
360#endif
361  UINT128 res;
362
363  UINT128 sign;
364  UINT32 comb, exp;
365  UINT128 trailing;
366  UINT128 d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11;
367  UINT128 bcoeff;
368  UINT128 dcoeff;
369  UINT64 nanb = 0;
370
371  sign.w[1] = (ba.w[HIGH_128W] & 0x8000000000000000ull);
372  sign.w[0] = 0;
373  comb = (ba.w[HIGH_128W] & 0x7fffc00000000000ull) >> 46;
374  trailing.w[1] = (ba.w[HIGH_128W] & 0x00003fffffffffffull);
375  trailing.w[0] = ba.w[LOW_128W];
376  exp = 0;
377
378  if ((comb & 0x1f000) == 0x1e000) {	// G0..G4 = 1111 -> Inf
379    res.w[HIGH_128W] = ba.w[HIGH_128W] & 0xf800000000000000ull;
380    res.w[LOW_128W] = 0;
381    BID_RETURN (res);
382    // Detect NaN, and canonicalize trailing
383  } else if ((comb & 0x1f000) == 0x1f000) {
384    if ((trailing.w[1] > 0x0000314dc6448d93ULL) ||	// significand is non-canonical
385	((trailing.w[1] == 0x0000314dc6448d93ULL)
386	 && (trailing.w[0] >= 0x38c15b0a00000000ULL))
387	// significand is non-canonical
388      ) {
389      trailing.w[1] = trailing.w[0] = 0ull;
390    }
391    bcoeff.w[1] = trailing.w[1];
392    bcoeff.w[0] = trailing.w[0];
393    nanb = ba.w[HIGH_128W] & 0xfe00000000000000ull;
394    exp = 0;
395  } else {	// Normal number
396    if ((comb & 0x18000) == 0x18000) {	// G0..G1 = 11 -> exp is G2..G11
397      exp = (comb >> 1) & 0x3fff;
398      bcoeff.w[1] =
399	((UINT64) (8 + (comb & 1)) << (UINT64) 46) | trailing.w[1];
400      bcoeff.w[0] = trailing.w[0];
401    } else {
402      exp = (comb >> 3) & 0x3fff;
403      bcoeff.w[1] =
404	((UINT64) (comb & 7) << (UINT64) 46) | trailing.w[1];
405      bcoeff.w[0] = trailing.w[0];
406    }
407    // Zero the coefficient if non-canonical (>= 10^34)
408    if (bcoeff.w[1] > 0x1ed09bead87c0ull ||
409	(bcoeff.w[1] == 0x1ed09bead87c0ull
410	 && bcoeff.w[0] >= 0x378D8E6400000000ull)) {
411      bcoeff.w[0] = bcoeff.w[1] = 0;
412    }
413  }
414  // Constant 2^128 / 1000 + 1
415  {
416    UINT128 t;
417    UINT64 t2;
418    UINT128 d1000;
419    UINT128 b11, b10, b9, b8, b7, b6, b5, b4, b3, b2, b1;
420    d1000.w[1] = 0x4189374BC6A7EFull;
421    d1000.w[0] = 0x9DB22D0E56041894ull;
422    __mul_128x128_high (b11, bcoeff, d1000);
423    __mul_128x128_high (b10, b11, d1000);
424    __mul_128x128_high (b9, b10, d1000);
425    __mul_128x128_high (b8, b9, d1000);
426    __mul_128x128_high (b7, b8, d1000);
427    __mul_128x128_high (b6, b7, d1000);
428    __mul_128x128_high (b5, b6, d1000);
429    __mul_128x128_high (b4, b5, d1000);
430    __mul_128x128_high (b3, b4, d1000);
431    __mul_128x128_high (b2, b3, d1000);
432    __mul_128x128_high (b1, b2, d1000);
433
434
435    __mul_64x128_full (t2, t, 1000ull, b11);
436    __sub_128_128 (d11, bcoeff, t);
437    __mul_64x128_full (t2, t, 1000ull, b10);
438    __sub_128_128 (d10, b11, t);
439    __mul_64x128_full (t2, t, 1000ull, b9);
440    __sub_128_128 (d9, b10, t);
441    __mul_64x128_full (t2, t, 1000ull, b8);
442    __sub_128_128 (d8, b9, t);
443    __mul_64x128_full (t2, t, 1000ull, b7);
444    __sub_128_128 (d7, b8, t);
445    __mul_64x128_full (t2, t, 1000ull, b6);
446    __sub_128_128 (d6, b7, t);
447    __mul_64x128_full (t2, t, 1000ull, b5);
448    __sub_128_128 (d5, b6, t);
449    __mul_64x128_full (t2, t, 1000ull, b4);
450    __sub_128_128 (d4, b5, t);
451    __mul_64x128_full (t2, t, 1000ull, b3);
452    __sub_128_128 (d3, b4, t);
453    __mul_64x128_full (t2, t, 1000ull, b2);
454    __sub_128_128 (d2, b3, t);
455    __mul_64x128_full (t2, t, 1000ull, b1);
456    __sub_128_128 (d1, b2, t);
457    d0 = b1;
458
459  }
460
461  dcoeff.w[0] = b2d[d11.w[0]] | (b2d[d10.w[0]] << 10) |
462    (b2d[d9.w[0]] << 20) | (b2d[d8.w[0]] << 30) | (b2d[d7.w[0]] << 40) |
463    (b2d[d6.w[0]] << 50) | (b2d[d5.w[0]] << 60);
464  dcoeff.w[1] =
465    (b2d[d5.w[0]] >> 4) | (b2d[d4.w[0]] << 6) | (b2d[d3.w[0]] << 16) |
466    (b2d[d2.w[0]] << 26) | (b2d[d1.w[0]] << 36);
467
468  res.w[0] = dcoeff.w[0];
469  if (d0.w[0] >= 8) {
470    res.w[1] =
471      sign.
472      w[1] |
473      ((0x18000 | ((exp >> 12) << 13) | ((d0.w[0] & 1) << 12) |
474	(exp & 0xfff)) << 46) | dcoeff.w[1];
475  } else {
476    res.w[1] =
477      sign.
478      w[1] | ((((exp >> 12) << 15) | (d0.w[0] << 12) | (exp & 0xfff))
479	      << 46) | dcoeff.w[1];
480  }
481
482  res.w[1] |= nanb;
483
484  BID_SWAP128 (res);
485  BID_RETURN (res);
486}
487
488#if DECIMAL_CALL_BY_REFERENCE
489void
490dpd_to_bid128 (UINT128 * pres, UINT128 * pda) {
491  UINT128 da = *pda;
492#else
493UINT128
494dpd_to_bid128 (UINT128 da) {
495#endif
496  UINT128 in = *(UINT128 *) & da;
497  UINT128 res;
498
499  UINT128 sign;
500  UINT64 exp, comb;
501  UINT128 trailing;
502  UINT64 d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11;
503  UINT128 bcoeff;
504  UINT64 tl, th;
505  UINT64 nanb = 0;
506
507  sign.w[1] = (in.w[HIGH_128W] & 0x8000000000000000ull);
508  sign.w[0] = 0;
509  comb = (in.w[HIGH_128W] & 0x7fffc00000000000ull) >> 46;
510  trailing.w[1] = (in.w[HIGH_128W] & 0x00003fffffffffffull);
511  trailing.w[0] = in.w[LOW_128W];
512  exp = 0;
513
514  if ((comb & 0x1f000) == 0x1e000) {	// G0..G4 = 1111 -> Inf
515    res.w[HIGH_128W] = in.w[HIGH_128W] & 0xf800000000000000ull;
516    res.w[LOW_128W] = 0ull;
517    BID_RETURN (res);
518  } else if ((comb & 0x1f000) == 0x1f000) {	// G0..G4 = 11111 -> NaN
519    nanb = in.w[HIGH_128W] & 0xfe00000000000000ull;
520    exp = 0;
521    d0 = 0;
522  } else {	// Normal number
523    if ((comb & 0x18000) == 0x18000) {	// G0..G1 = 11 -> d0 = 8 + G4
524      d0 = 8 + (comb & 0x01000 ? 1 : 0);
525      exp =
526	(comb & 0x04000 ? 1 : 0) * 0x2000 +
527	(comb & 0x02000 ? 1 : 0) * 0x1000;
528      // exp leading bits are G2..G3
529    } else {
530      d0 =
531	4 * (comb & 0x04000 ? 1 : 0) + 2 * (comb & 0x2000 ? 1 : 0) +
532	(comb & 0x1000 ? 1 : 0);
533      exp =
534	(comb & 0x10000 ? 1 : 0) * 0x2000 +
535	(comb & 0x08000 ? 1 : 0) * 0x1000;
536      // exp loading bits are G0..G1
537    }
538  }
539
540  d11 = d2b[(trailing.w[0]) & 0x3ff];
541  d10 = d2b[(trailing.w[0] >> 10) & 0x3ff];
542  d9 = d2b[(trailing.w[0] >> 20) & 0x3ff];
543  d8 = d2b[(trailing.w[0] >> 30) & 0x3ff];
544  d7 = d2b[(trailing.w[0] >> 40) & 0x3ff];
545  d6 = d2b[(trailing.w[0] >> 50) & 0x3ff];
546  d5 = d2b[(trailing.w[0] >> 60) | ((trailing.w[1] & 0x3f) << 4)];
547  d4 = d2b[(trailing.w[1] >> 6) & 0x3ff];
548  d3 = d2b[(trailing.w[1] >> 16) & 0x3ff];
549  d2 = d2b[(trailing.w[1] >> 26) & 0x3ff];
550  d1 = d2b[(trailing.w[1] >> 36) & 0x3ff];
551
552  tl =
553    d11 + (d10 * 1000ull) + (d9 * 1000000ull) + (d8 * 1000000000ull) +
554    (d7 * 1000000000000ull) + (d6 * 1000000000000000ull);
555  th =
556    d5 + (d4 * 1000ull) + (d3 * 1000000ull) + (d2 * 1000000000ull) +
557    (d1 * 1000000000000ull) + (d0 * 1000000000000000ull);
558  __mul_64x64_to_128 (bcoeff, th, 1000000000000000000ull);
559  __add_128_64 (bcoeff, bcoeff, tl);
560
561  if (!nanb)
562    exp += (comb & 0xfff);
563
564  res.w[0] = bcoeff.w[0];
565  res.w[1] = (exp << 49) | sign.w[1] | bcoeff.w[1];
566
567  res.w[1] |= nanb;
568
569  BID_SWAP128 (res);
570  BID_RETURN (res);
571}
572
573UINT128
574bid_to_bid128 (UINT128 bq) {
575  UINT128 res;
576  UINT64 sign, comb, exp;
577  UINT64 trailing;
578  UINT64 bcoeff;
579
580  UINT128 rq;
581  UINT128 qcoeff;
582  UINT64 ba, bb;
583
584  ba = *((UINT64 *) & bq + HIGH_128W);
585  bb = *((UINT64 *) & bq + LOW_128W);
586
587  sign = (ba & 0x8000000000000000ull);
588  comb = (ba & 0x7fffc00000000000ull) >> 46;
589  trailing = (ba & 0x00003fffffffffffull);
590
591  if ((comb & 0x18000) == 0x18000) {	// G0..G1 = 11 -> exp is G2..G11
592    exp = (comb >> 1) & 0x3fff;
593    bcoeff = ((8 + (comb & 1)) << 46) | trailing;
594  } else {
595    exp = (comb >> 3) & 0x3fff;
596    bcoeff = ((comb & 7) << 46) | trailing;
597  }
598
599  if ((comb & 0x1f000) == 0x1f000) {	//NaN
600    ba &= 0xfe003fffffffffffULL;	// make exponent 0
601    bcoeff &= 0x00003fffffffffffull;	// NaN payloat is only T.
602    if ((bcoeff > 0x0000314dc6448d93ULL) ||	// significand is non-canonical
603	((bcoeff == 0x0000314dc6448d93ULL)
604	 && (bb >= 0x38c15b0a00000000ULL))
605	// significand is non-canonical
606      ) {
607      bcoeff = 0ull;
608      ba &= ~0x00003fffffffffffull;
609      bb = 0ull;
610    }
611    *((UINT64 *) & rq + HIGH_128W) = ba;
612    *((UINT64 *) & rq + LOW_128W) = bb;
613    return rq;
614  } else if ((comb & 0x1e000) == 0x1e000) {	//Inf
615    ba &= 0xf800000000000000ULL;	// make exponent and significand 0
616    bb = 0;
617    *((UINT64 *) & rq + HIGH_128W) = ba;
618    *((UINT64 *) & rq + LOW_128W) = bb;
619    return rq;
620  }
621
622  if ((bcoeff > 0x0001ed09bead87c0ull)
623      || ((bcoeff == 0x0001ed09bead87c0ull)
624	  && (bb > 0x378d8e63ffffffffull))) {
625    // significand is non-canonical
626    bcoeff = 0ull;
627    bb = 0ull;
628  }
629
630  *((UINT64 *) & qcoeff + 1) = bcoeff;
631  *((UINT64 *) & qcoeff + 0) = bb;
632
633  get_BID128_fast (&res, sign, exp, qcoeff);
634
635  BID_SWAP128 (res);
636  return res;
637}
638
639UINT32
640bid32_canonize (UINT32 ba) {
641  FPSC bidrnd;
642  unsigned int rnd = 0;
643
644  UINT32 res;
645  UINT32 sign, comb, exp;
646  UINT32 trailing;
647  UINT32 bcoeff;
648
649  sign = (ba & 0x80000000);
650  comb = (ba & 0x7ff00000) >> 20;
651  trailing = (ba & 0x000fffff);
652
653  if ((comb & 0x600) == 0x600) {	// G0..G1 = 11 -> exp is G2..G11
654    exp = (comb >> 1) & 0xff;
655    bcoeff = ((8 + (comb & 1)) << 20) | trailing;
656  } else {
657    exp = (comb >> 3) & 0xff;
658    bcoeff = ((comb & 7) << 20) | trailing;
659  }
660
661  if ((comb & 0x7c0) == 0x7c0) {	//NaN
662    ba &= 0xfe0fffff;	// make exponent 0
663    bcoeff &= 0x000fffff;	// NaN payloat is only T.
664    if (bcoeff >= 1000000)
665      ba &= 0xfff00000;	//treat non-canonical significand
666    return ba;
667  } else if ((comb & 0x780) == 0x780) {	//Inf
668    ba &= 0xf8000000;	// make exponent and significand 0
669    return ba;
670  }
671
672  if (bcoeff >= 10000000)
673    bcoeff = 0;
674  rnd = bidrnd = ROUNDING_TO_NEAREST;
675  res = get_BID32 (sign, exp, bcoeff, rnd, &bidrnd);
676  return res;
677}
678
679UINT64
680bid64_canonize (UINT64 ba) {
681  UINT64 res;
682  UINT64 sign, comb, exp;
683  UINT64 trailing;
684  UINT64 bcoeff;
685
686  sign = (ba & 0x8000000000000000ull);
687  comb = (ba & 0x7ffc000000000000ull) >> 50;
688  trailing = (ba & 0x0003ffffffffffffull);
689
690
691  if ((comb & 0x1800) == 0x1800) {	// G0..G1 = 11 -> exp is G2..G11
692    exp = (comb >> 1) & 0x3ff;
693    bcoeff = ((8 + (comb & 1)) << 50) | trailing;
694  } else {
695    exp = (comb >> 3) & 0x3ff;
696    bcoeff = ((comb & 7) << 50) | trailing;
697  }
698
699  if ((comb & 0x1f00) == 0x1f00) {	//NaN
700    ba &= 0xfe03ffffffffffffULL;	// make exponent 0
701    bcoeff &= 0x0003ffffffffffffull;	// NaN payloat is only T.
702    if (bcoeff >= 1000000000000000ull)
703      ba &= 0xfe00000000000000ull;	// treat non canonical significand and zero G6-G12
704    return ba;
705  } else if ((comb & 0x1e00) == 0x1e00) {	//Inf
706    ba &= 0xf800000000000000ULL;	// make exponent and significand 0
707    return ba;
708  }
709
710  if (bcoeff >= 10000000000000000ull) {
711    bcoeff = 0ull;
712  }
713  res = very_fast_get_BID64 (sign, exp, bcoeff);
714  return res;
715}
716
717UINT128
718bid128_canonize (UINT128 bq) {
719  UINT128 res;
720  UINT64 sign, comb, exp;
721  UINT64 trailing;
722  UINT64 bcoeff;
723
724  UINT128 rq;
725  UINT128 qcoeff;
726  UINT64 ba, bb;
727
728  ba = *((UINT64 *) & bq + HIGH_128W);
729  bb = *((UINT64 *) & bq + LOW_128W);
730
731  sign = (ba & 0x8000000000000000ull);
732  comb = (ba & 0x7fffc00000000000ull) >> 46;
733  trailing = (ba & 0x00003fffffffffffull);
734
735  if ((comb & 0x18000) == 0x18000) {	// G0..G1 = 11 -> exp is G2..G11
736    exp = (comb >> 1) & 0x3fff;
737    bcoeff = ((8 + (comb & 1)) << 46) | trailing;
738  } else {
739    exp = (comb >> 3) & 0x3fff;
740    bcoeff = ((comb & 7) << 46) | trailing;
741  }
742
743  if ((comb & 0x1f000) == 0x1f000) {	//NaN
744    ba &= 0xfe003fffffffffffULL;	// make exponent 0
745    bcoeff &= 0x00003fffffffffffull;	// NaN payload is only T.
746
747    if ((bcoeff > 0x0000314dc6448d93ULL) ||	// significand is non-canonical
748	((bcoeff == 0x0000314dc6448d93ULL)
749	 && (bb >= 0x38c15b0a00000000ULL))
750	// significand is non-canonical
751      ) {
752      bcoeff = 0ull;
753      ba &= ~0x00003fffffffffffull;
754      bb = 0ull;
755    }
756    *((UINT64 *) & rq + HIGH_128W) = ba;
757    *((UINT64 *) & rq + LOW_128W) = bb;
758    return rq;
759  } else if ((comb & 0x1e000) == 0x1e000) {	//Inf
760    ba &= 0xf800000000000000ULL;	// make exponent and significand 0
761    bb = 0;
762    *((UINT64 *) & rq + HIGH_128W) = ba;
763    *((UINT64 *) & rq + LOW_128W) = bb;
764    return rq;
765  }
766
767  if ((bcoeff > 0x0001ed09bead87c0ull) ||	// significand is non-canonical
768      ((bcoeff == 0x0001ed09bead87c0ull)
769       && (bb > 0x378d8e63ffffffffull))
770      // significand is non-canonical
771    ) {
772    bcoeff = 0ull;
773    bb = 0ull;
774  }
775
776  *((UINT64 *) & qcoeff + 1) = bcoeff;
777  *((UINT64 *) & qcoeff + 0) = bb;
778
779  get_BID128_fast (&res, sign, exp, qcoeff);
780  BID_SWAP128 (res);
781  return res;
782}
783