1/* Copyright (C) 2007-2022 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#undef IN_LIBGCC2
25#include "bid-dpd.h"
26
27/* get full 64x64bit product */
28#define __mul_64x64_to_128(P, CX, CY)             \
29{                                                 \
30  UINT64 CXH, CXL, CYH,CYL,PL,PH,PM,PM2;  \
31  CXH = (CX) >> 32;                               \
32  CXL = (UINT32)(CX);                             \
33  CYH = (CY) >> 32;                               \
34  CYL = (UINT32)(CY);                             \
35                                                  \
36  PM = CXH*CYL;                                   \
37  PH = CXH*CYH;                                   \
38  PL = CXL*CYL;                                   \
39  PM2 = CXL*CYH;                                  \
40  PH += (PM>>32);                                 \
41  PM = (UINT64)((UINT32)PM)+PM2+(PL>>32);         \
42                                                  \
43  (P).w[1] = PH + (PM>>32);                       \
44  (P).w[0] = (PM<<32)+(UINT32)PL;                 \
45}
46
47/* add 64-bit value to 128-bit */
48#define __add_128_64(R128, A128, B64)             \
49{                                                 \
50  UINT64 R64H;                                    \
51  R64H = (A128).w[1];                             \
52  (R128).w[0] = (B64) + (A128).w[0];              \
53  if((R128).w[0] < (B64)) R64H ++;                \
54  (R128).w[1] = R64H;                             \
55}
56
57/* add 128-bit value to 128-bit (assume no carry-out) */
58#define __add_128_128(R128, A128, B128)           \
59{                                                 \
60  UINT128 Q128;                                   \
61  Q128.w[1] = (A128).w[1]+(B128).w[1];            \
62  Q128.w[0] = (B128).w[0] + (A128).w[0];          \
63  if(Q128.w[0] < (B128).w[0]) Q128.w[1] ++;       \
64  (R128).w[1] = Q128.w[1];                        \
65  (R128).w[0] = Q128.w[0];                        \
66}
67
68#define __mul_128x128_high(Q, A, B)               \
69{                                                 \
70  UINT128 ALBL, ALBH, AHBL, AHBH, QM, QM2;        \
71                                                  \
72  __mul_64x64_to_128(ALBH, (A).w[0], (B).w[1]);   \
73  __mul_64x64_to_128(AHBL, (B).w[0], (A).w[1]);   \
74  __mul_64x64_to_128(ALBL, (A).w[0], (B).w[0]);   \
75  __mul_64x64_to_128(AHBH, (A).w[1],(B).w[1]);    \
76                                                  \
77  __add_128_128(QM, ALBH, AHBL);                  \
78  __add_128_64(QM2, QM, ALBL.w[1]);               \
79  __add_128_64((Q), AHBH, QM2.w[1]);              \
80}
81
82#include "bid2dpd_dpd2bid.h"
83
84static const unsigned int dm103[] =
85  { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 11000 };
86
87void _bid_to_dpd32 (_Decimal32 *, _Decimal32 *);
88
89void
90_bid_to_dpd32 (_Decimal32 *pres, _Decimal32 *px) {
91  unsigned int sign, coefficient_x, exp, dcoeff;
92  unsigned int b2, b1, b0, b01, res;
93  _Decimal32 x = *px;
94
95  sign = (x & 0x80000000);
96  if ((x & 0x60000000ul) == 0x60000000ul) {
97    /* special encodings */
98    if ((x & 0x78000000ul) == 0x78000000ul) {
99      *pres = x; /* NaN or Infinity */
100      return;
101    }
102    /* coefficient */
103    coefficient_x = (x & 0x001ffffful) | 0x00800000ul;
104    if (coefficient_x >= 10000000) coefficient_x = 0;
105    /* get exponent */
106    exp = (x >> 21) & 0xff;
107  } else {
108    exp = (x >> 23) & 0xff;
109    coefficient_x = (x & 0x007ffffful);
110  }
111  b01 = coefficient_x / 1000;
112  b2 = coefficient_x - 1000 * b01;
113  b0 = b01 / 1000;
114  b1 = b01 - 1000 * b0;
115  dcoeff = b2d[b2] | b2d2[b1];
116  if (b0 >= 8) { /* is b0 8 or 9? */
117    res = sign | ((0x600 | ((exp >> 6) << 7) |
118        ((b0 & 1) << 6) | (exp & 0x3f)) << 20) | dcoeff;
119  } else { /* else b0 is 0..7 */
120    res = sign | ((((exp >> 6) << 9) | (b0 << 6) |
121        (exp & 0x3f)) << 20) | dcoeff;
122  }
123  *pres = res;
124}
125
126void _dpd_to_bid32 (_Decimal32 *, _Decimal32 *);
127
128void
129_dpd_to_bid32 (_Decimal32 *pres, _Decimal32 *px) {
130  unsigned int r;
131  unsigned int sign, exp, bcoeff;
132  UINT64 trailing;
133  unsigned int d0, d1, d2;
134  _Decimal32 x = *px;
135
136  sign = (x & 0x80000000);
137  trailing = (x & 0x000fffff);
138  if ((x & 0x78000000) == 0x78000000) {
139    *pres = x;
140    return;
141  }
142  /* normal number */
143  if ((x & 0x60000000) == 0x60000000) { /* G0..G1 = 11 -> d0 = 8 + G4 */
144    d0 = d2b3[((x >> 26) & 1) | 8]; /* d0 = (comb & 0x0100 ? 9 : 8); */
145    exp = (x >> 27) & 3; /* exp leading bits are G2..G3 */
146  } else {
147    d0 = d2b3[(x >> 26) & 0x7];
148    exp = (x >> 29) & 3; /* exp loading bits are G0..G1 */
149  }
150  d1 = d2b2[(trailing >> 10) & 0x3ff];
151  d2 = d2b[(trailing) & 0x3ff];
152  bcoeff = d2 + d1 + d0;
153  exp = (exp << 6) + ((x >> 20) & 0x3f);
154  if (bcoeff < (1 << 23)) {
155    r = exp;
156    r <<= 23;
157    r |= (bcoeff | sign);
158  } else {
159    r = exp;
160    r <<= 21;
161    r |= (sign | 0x60000000ul);
162    /* add coeff, without leading bits */
163    r |= (((unsigned int) bcoeff) & 0x1fffff);
164  }
165  *pres = r;
166}
167
168void _bid_to_dpd64 (_Decimal64 *, _Decimal64 *);
169
170void
171_bid_to_dpd64 (_Decimal64 *pres, _Decimal64 *px) {
172  UINT64 res;
173  UINT64 sign, comb, exp, B34, B01;
174  UINT64 d103, D61;
175  UINT64 b0, b2, b3, b5;
176  unsigned int b1, b4;
177  UINT64 bcoeff;
178  UINT64 dcoeff;
179  unsigned int yhi, ylo;
180  _Decimal64 x = *px;
181
182  sign = (x & 0x8000000000000000ull);
183  comb = (x & 0x7ffc000000000000ull) >> 51;
184  if ((comb & 0xf00) == 0xf00) {
185    *pres = x;
186    return;
187  }
188  /* Normal number */
189  if ((comb & 0xc00) == 0xc00) { /* G0..G1 = 11 -> exp is G2..G11 */
190    exp = (comb) & 0x3ff;
191    bcoeff = (x & 0x0007ffffffffffffull) | 0x0020000000000000ull;
192    if (bcoeff >= 10000000000000000ull)
193      bcoeff = 0;
194  } else {
195    exp = (comb >> 2) & 0x3ff;
196    bcoeff = (x & 0x001fffffffffffffull);
197  }
198  D61 = 2305843009ull; /* Floor(2^61 / 10^9) */
199  /* Multiply the binary coefficient by ceil(2^64 / 1000), and take the upper
200     64-bits in order to compute a division by 1000. */
201  yhi = (D61 * (UINT64)(bcoeff >> (UINT64)27)) >> (UINT64)34;
202  ylo = bcoeff - 1000000000ull * yhi;
203  if (ylo >= 1000000000) {
204    ylo = ylo - 1000000000;
205    yhi = yhi + 1;
206  }
207  d103 = 0x4189374c;
208  B34 = ((UINT64) ylo * d103) >> (32 + 8);
209  B01 = ((UINT64) yhi * d103) >> (32 + 8);
210  b5 = ylo - B34 * 1000;
211  b2 = yhi - B01 * 1000;
212  b3 = ((UINT64) B34 * d103) >> (32 + 8);
213  b0 = ((UINT64) B01 * d103) >> (32 + 8);
214  b4 = (unsigned int) B34 - (unsigned int) b3 *1000;
215  b1 = (unsigned int) B01 - (unsigned int) dm103[b0];
216  dcoeff = b2d[b5] | b2d2[b4] | b2d3[b3] | b2d4[b2] | b2d5[b1];
217  if (b0 >= 8) /* is b0 8 or 9? */
218    res = sign | ((0x1800 | ((exp >> 8) << 9) | ((b0 & 1) << 8) |
219                   (exp & 0xff)) << 50) | dcoeff;
220  else /* else b0 is 0..7 */
221    res = sign | ((((exp >> 8) << 11) | (b0 << 8) |
222                     (exp & 0xff)) << 50) | dcoeff;
223  *pres = res;
224}
225
226void _dpd_to_bid64 (_Decimal64 *, _Decimal64 *);
227
228void
229_dpd_to_bid64 (_Decimal64 *pres, _Decimal64 *px) {
230  UINT64 res;
231  UINT64 sign, comb, exp;
232  UINT64 trailing;
233  UINT64 d0, d1, d2;
234  unsigned int d3, d4, d5;
235  UINT64 bcoeff, mask;
236  _Decimal64 x = *px;
237
238  sign = (x & 0x8000000000000000ull);
239  comb = (x & 0x7ffc000000000000ull) >> 50;
240  trailing = (x & 0x0003ffffffffffffull);
241  if ((comb & 0x1e00) == 0x1e00) {
242    *pres = x;
243    return;
244  }
245  /* normal number */
246  if ((comb & 0x1800) == 0x1800) { /* G0..G1 = 11 -> d0 = 8 + G4 */
247    d0 = d2b6[((comb >> 8) & 1) | 8]; /* d0 = (comb & 0x0100 ? 9 : 8); */
248    exp = (comb & 0x600) >> 1; /* exp = (comb & 0x0400 ? 1 : 0) * 0x200 +
249        (comb & 0x0200 ? 1 : 0) * 0x100; exp leading bits are G2..G3 */
250  } else {
251    d0 = d2b6[(comb >> 8) & 0x7];
252    exp = (comb & 0x1800) >> 3; /* exp = (comb & 0x1000 ? 1 : 0) * 0x200 +
253        (comb & 0x0800 ? 1 : 0) * 0x100; exp loading bits are G0..G1 */
254  }
255  d1 = d2b5[(trailing >> 40) & 0x3ff];
256  d2 = d2b4[(trailing >> 30) & 0x3ff];
257  d3 = d2b3[(trailing >> 20) & 0x3ff];
258  d4 = d2b2[(trailing >> 10) & 0x3ff];
259  d5 = d2b[(trailing) & 0x3ff];
260  bcoeff = (d5 + d4 + d3) + d2 + d1 + d0;
261  exp += (comb & 0xff);
262  mask = 1;
263  mask <<= 53;
264  if (bcoeff < mask) { /* check whether coefficient fits in 10*5+3 bits */
265    res = exp;
266    res <<= 53;
267    res |= (bcoeff | sign);
268    *pres = res;
269    return;
270  }
271  /* special format */
272  res = (exp << 51) | (sign | 0x6000000000000000ull);
273  /* add coeff, without leading bits */
274  mask = (mask >> 2) - 1;
275  bcoeff &= mask;
276  res |= bcoeff;
277  *pres = res;
278}
279
280void _bid_to_dpd128 (_Decimal128 *, _Decimal128 *);
281
282void
283_bid_to_dpd128 (_Decimal128 *pres, _Decimal128 *px) {
284  UINT128 res;
285  UINT128 sign;
286  unsigned int comb;
287  UINT128 bcoeff;
288  UINT128 dcoeff;
289  UINT128 BH, d1018, BT2, BT1;
290  UINT64 exp, BL, d109;
291  UINT64 d106, d103;
292  UINT64 k1, k2, k4, k5, k7, k8, k10, k11;
293  unsigned int BHH32, BLL32, BHL32, BLH32, k0, k3, k6, k9, amount;
294  _Decimal128 x = *px;
295
296  sign.w[1] = (x.w[1] & 0x8000000000000000ull);
297  sign.w[0] = 0;
298  comb = (x.w[1] /*& 0x7fffc00000000000ull */ ) >> 46;
299  exp = 0;
300  if ((comb & 0x1e000) == 0x1e000) {
301    res = x;
302  } else { /* normal number */
303    if ((comb & 0x18000) == 0x18000) {
304      /* Noncanonical significand (prepending 8 or 9 to any 110-bit
305	 trailing significand field produces a value above 10^34).  */
306      exp = (comb & 0x7fff) >> 1;
307      bcoeff.w[1] = 0;
308      bcoeff.w[0] = 0;
309    } else {
310      exp = ((x.w[1] & 0x7fff000000000000ull) >> 49) & 0x3fff;
311      bcoeff.w[1] = (x.w[1] & 0x0001ffffffffffffull);
312      bcoeff.w[0] = x.w[0];
313      if (bcoeff.w[1] > 0x1ed09bead87c0ull
314	  || (bcoeff.w[1] == 0x1ed09bead87c0ull
315	      && bcoeff.w[0] >= 0x378d8e6400000000ull)) {
316	bcoeff.w[1] = 0;
317	bcoeff.w[0] = 0;
318      }
319    }
320    d1018 = reciprocals10_128[18];
321    __mul_128x128_high (BH, bcoeff, d1018);
322    amount = recip_scale[18];
323    BH.w[0] = (BH.w[0] >> amount) | (BH.w[1] << (64 - amount));
324    BL = bcoeff.w[0] - BH.w[0] * 1000000000000000000ull;
325    d109 = reciprocals10_64[9];
326    __mul_64x64_to_128 (BT1, BH.w[0], d109);
327    BHH32 = (unsigned int) (BT1.w[1] >> short_recip_scale[9]);
328    BHL32 = (unsigned int) BH.w[0] - BHH32 * 1000000000;
329    __mul_64x64_to_128 (BT2, BL, d109);
330    BLH32 = (unsigned int) (BT2.w[1] >> short_recip_scale[9]);
331    BLL32 = (unsigned int) BL - BLH32 * 1000000000;
332    d106 = 0x431BDE83;
333    d103 = 0x4189374c;
334    k0 = ((UINT64) BHH32 * d106) >> (32 + 18);
335    BHH32 -= (unsigned int) k0 *1000000;
336    k1 = ((UINT64) BHH32 * d103) >> (32 + 8);
337    k2 = BHH32 - (unsigned int) k1 *1000;
338    k3 = ((UINT64) BHL32 * d106) >> (32 + 18);
339    BHL32 -= (unsigned int) k3 *1000000;
340    k4 = ((UINT64) BHL32 * d103) >> (32 + 8);
341    k5 = BHL32 - (unsigned int) k4 *1000;
342    k6 = ((UINT64) BLH32 * d106) >> (32 + 18);
343    BLH32 -= (unsigned int) k6 *1000000;
344    k7 = ((UINT64) BLH32 * d103) >> (32 + 8);
345    k8 = BLH32 - (unsigned int) k7 *1000;
346    k9 = ((UINT64) BLL32 * d106) >> (32 + 18);
347    BLL32 -= (unsigned int) k9 *1000000;
348    k10 = ((UINT64) BLL32 * d103) >> (32 + 8);
349    k11 = BLL32 - (unsigned int) k10 *1000;
350    dcoeff.w[1] = (b2d[k5] >> 4) | (b2d[k4] << 6) | (b2d[k3] << 16) |
351      (b2d[k2] << 26) | (b2d[k1] << 36);
352    dcoeff.w[0] = b2d[k11] | (b2d[k10] << 10) | (b2d[k9] << 20) |
353      (b2d[k8] << 30) | (b2d[k7] << 40) | (b2d[k6] << 50) | (b2d[k5] << 60);
354    res.w[0] = dcoeff.w[0];
355    if (k0 >= 8) {
356      res.w[1] = sign.w[1] | ((0x18000 | ((exp >> 12) << 13) |
357          ((k0 & 1) << 12) | (exp & 0xfff)) << 46) | dcoeff.w[1];
358    } else {
359      res.w[1] = sign.w[1] | ((((exp >> 12) << 15) | (k0 << 12) |
360          (exp & 0xfff)) << 46) | dcoeff.w[1];
361    }
362  }
363  *pres = res;
364}
365
366void _dpd_to_bid128 (_Decimal128 *, _Decimal128 *);
367
368void
369_dpd_to_bid128 (_Decimal128 *pres, _Decimal128 *px) {
370  UINT128 res;
371  UINT128 sign;
372  UINT64 exp, comb;
373  UINT128 trailing;
374  UINT64 d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11;
375  UINT128 bcoeff;
376  UINT64 tl, th;
377  _Decimal128 x = *px;
378
379  sign.w[1] = (x.w[1] & 0x8000000000000000ull);
380  sign.w[0] = 0;
381  comb = (x.w[1] & 0x7fffc00000000000ull) >> 46;
382  trailing.w[1] = x.w[1];
383  trailing.w[0] = x.w[0];
384  if ((comb & 0x1e000) == 0x1e000) {
385      *pres = x;
386      return;
387  }
388  if ((comb & 0x18000) == 0x18000) { /* G0..G1 = 11 -> d0 = 8 + G4 */
389    d0 = d2b6[8 + ((comb & 0x01000) >> 12)];
390    exp = (comb & 0x06000) >> 1;  /* exp leading bits are G2..G3 */
391  } else {
392    d0 = d2b6[((comb & 0x07000) >> 12)];
393    exp = (comb & 0x18000) >> 3;  /* exp loading bits are G0..G1 */
394  }
395  d11 = d2b[(trailing.w[0]) & 0x3ff];
396  d10 = d2b2[(trailing.w[0] >> 10) & 0x3ff];
397  d9 = d2b3[(trailing.w[0] >> 20) & 0x3ff];
398  d8 = d2b4[(trailing.w[0] >> 30) & 0x3ff];
399  d7 = d2b5[(trailing.w[0] >> 40) & 0x3ff];
400  d6 = d2b6[(trailing.w[0] >> 50) & 0x3ff];
401  d5 = d2b[(trailing.w[0] >> 60) | ((trailing.w[1] & 0x3f) << 4)];
402  d4 = d2b2[(trailing.w[1] >> 6) & 0x3ff];
403  d3 = d2b3[(trailing.w[1] >> 16) & 0x3ff];
404  d2 = d2b4[(trailing.w[1] >> 26) & 0x3ff];
405  d1 = d2b5[(trailing.w[1] >> 36) & 0x3ff];
406  tl = d11 + d10 + d9 + d8 + d7 + d6;
407  th = d5 + d4 + d3 + d2 + d1 + d0;
408  __mul_64x64_to_128 (bcoeff, th, 1000000000000000000ull);
409  __add_128_64 (bcoeff, bcoeff, tl);
410  exp += (comb & 0xfff);
411  res.w[0] = bcoeff.w[0];
412  res.w[1] = (exp << 49) | sign.w[1] | bcoeff.w[1];
413  *pres = res;
414}
415