1// Copyright 2018 Ulf Adams
2//
3// The contents of this file may be used under the terms of the Apache License,
4// Version 2.0.
5//
6//    (See accompanying file LICENSE-Apache or copy at
7//     http://www.apache.org/licenses/LICENSE-2.0)
8//
9// Alternatively, the contents of this file may be used under the terms of
10// the Boost Software License, Version 1.0.
11//    (See accompanying file LICENSE-Boost or copy at
12//     https://www.boost.org/LICENSE_1_0.txt)
13//
14// Unless required by applicable law or agreed to in writing, this software
15// is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
16// KIND, either express or implied.
17#ifndef RYU_D2S_INTRINSICS_H
18#define RYU_D2S_INTRINSICS_H
19
20
21// Defines RYU_32_BIT_PLATFORM if applicable.
22
23// ABSL avoids uint128_t on Win32 even if __SIZEOF_INT128__ is defined.
24// Let's do the same for now.
25#if defined(__SIZEOF_INT128__) && !defined(_MSC_VER) && !defined(RYU_ONLY_64_BIT_OPS)
26#define HAS_UINT128
27#elif defined(_MSC_VER) && !defined(RYU_ONLY_64_BIT_OPS) && defined(_M_X64)
28#define HAS_64_BIT_INTRINSICS
29#endif
30
31#if defined(HAS_64_BIT_INTRINSICS)
32
33
34static inline uint64_t umul128(const uint64_t a, const uint64_t b, uint64_t* const productHi) {
35  return _umul128(a, b, productHi);
36}
37
38// Returns the lower 64 bits of (hi*2^64 + lo) >> dist, with 0 < dist < 64.
39static inline uint64_t shiftright128(const uint64_t lo, const uint64_t hi, const uint32_t dist) {
40  // For the __shiftright128 intrinsic, the shift value is always
41  // modulo 64.
42  // In the current implementation of the double-precision version
43  // of Ryu, the shift value is always < 64. (In the case
44  // RYU_OPTIMIZE_SIZE == 0, the shift value is in the range [49, 58].
45  // Otherwise in the range [2, 59].)
46  // However, this function is now also called by s2d, which requires supporting
47  // the larger shift range (TODO: what is the actual range?).
48  // Check this here in case a future change requires larger shift
49  // values. In this case this function needs to be adjusted.
50  assert(dist < 64);
51  return __shiftright128(lo, hi, (unsigned char) dist);
52}
53
54#else // defined(HAS_64_BIT_INTRINSICS)
55
56static inline uint64_t umul128(const uint64_t a, const uint64_t b, uint64_t* const productHi) {
57  // The casts here help MSVC to avoid calls to the __allmul library function.
58  const uint32_t aLo = (uint32_t)a;
59  const uint32_t aHi = (uint32_t)(a >> 32);
60  const uint32_t bLo = (uint32_t)b;
61  const uint32_t bHi = (uint32_t)(b >> 32);
62
63  const uint64_t b00 = (uint64_t)aLo * bLo;
64  const uint64_t b01 = (uint64_t)aLo * bHi;
65  const uint64_t b10 = (uint64_t)aHi * bLo;
66  const uint64_t b11 = (uint64_t)aHi * bHi;
67
68  const uint32_t b00Lo = (uint32_t)b00;
69  const uint32_t b00Hi = (uint32_t)(b00 >> 32);
70
71  const uint64_t mid1 = b10 + b00Hi;
72  const uint32_t mid1Lo = (uint32_t)(mid1);
73  const uint32_t mid1Hi = (uint32_t)(mid1 >> 32);
74
75  const uint64_t mid2 = b01 + mid1Lo;
76  const uint32_t mid2Lo = (uint32_t)(mid2);
77  const uint32_t mid2Hi = (uint32_t)(mid2 >> 32);
78
79  const uint64_t pHi = b11 + mid1Hi + mid2Hi;
80  const uint64_t pLo = ((uint64_t)mid2Lo << 32) | b00Lo;
81
82  *productHi = pHi;
83  return pLo;
84}
85
86static inline uint64_t shiftright128(const uint64_t lo, const uint64_t hi, const uint32_t dist) {
87  // We don't need to handle the case dist >= 64 here (see above).
88  assert(dist < 64);
89  assert(dist > 0);
90  return (hi << (64 - dist)) | (lo >> dist);
91}
92
93#endif // defined(HAS_64_BIT_INTRINSICS)
94
95#if defined(RYU_32_BIT_PLATFORM)
96
97// Returns the high 64 bits of the 128-bit product of a and b.
98static inline uint64_t umulh(const uint64_t a, const uint64_t b) {
99  // Reuse the umul128 implementation.
100  // Optimizers will likely eliminate the instructions used to compute the
101  // low part of the product.
102  uint64_t hi;
103  umul128(a, b, &hi);
104  return hi;
105}
106
107// On 32-bit platforms, compilers typically generate calls to library
108// functions for 64-bit divisions, even if the divisor is a constant.
109//
110// E.g.:
111// https://bugs.llvm.org/show_bug.cgi?id=37932
112// https://gcc.gnu.org/bugzilla/show_bug.cgi?id=17958
113// https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37443
114//
115// The functions here perform division-by-constant using multiplications
116// in the same way as 64-bit compilers would do.
117//
118// NB:
119// The multipliers and shift values are the ones generated by clang x64
120// for expressions like x/5, x/10, etc.
121
122static inline uint64_t div5(const uint64_t x) {
123  return umulh(x, 0xCCCCCCCCCCCCCCCDu) >> 2;
124}
125
126static inline uint64_t div10(const uint64_t x) {
127  return umulh(x, 0xCCCCCCCCCCCCCCCDu) >> 3;
128}
129
130static inline uint64_t div100(const uint64_t x) {
131  return umulh(x >> 2, 0x28F5C28F5C28F5C3u) >> 2;
132}
133
134static inline uint64_t div1e8(const uint64_t x) {
135  return umulh(x, 0xABCC77118461CEFDu) >> 26;
136}
137
138static inline uint64_t div1e9(const uint64_t x) {
139  return umulh(x >> 9, 0x44B82FA09B5A53u) >> 11;
140}
141
142static inline uint32_t mod1e9(const uint64_t x) {
143  // Avoid 64-bit math as much as possible.
144  // Returning (uint32_t) (x - 1000000000 * div1e9(x)) would
145  // perform 32x64-bit multiplication and 64-bit subtraction.
146  // x and 1000000000 * div1e9(x) are guaranteed to differ by
147  // less than 10^9, so their highest 32 bits must be identical,
148  // so we can truncate both sides to uint32_t before subtracting.
149  // We can also simplify (uint32_t) (1000000000 * div1e9(x)).
150  // We can truncate before multiplying instead of after, as multiplying
151  // the highest 32 bits of div1e9(x) can't affect the lowest 32 bits.
152  return ((uint32_t) x) - 1000000000 * ((uint32_t) div1e9(x));
153}
154
155#else // defined(RYU_32_BIT_PLATFORM)
156
157static inline uint64_t div5(const uint64_t x) {
158  return x / 5;
159}
160
161static inline uint64_t div10(const uint64_t x) {
162  return x / 10;
163}
164
165static inline uint64_t div100(const uint64_t x) {
166  return x / 100;
167}
168
169static inline uint64_t div1e8(const uint64_t x) {
170  return x / 100000000;
171}
172
173static inline uint64_t div1e9(const uint64_t x) {
174  return x / 1000000000;
175}
176
177static inline uint32_t mod1e9(const uint64_t x) {
178  return (uint32_t) (x - 1000000000 * div1e9(x));
179}
180
181#endif // defined(RYU_32_BIT_PLATFORM)
182
183static inline uint32_t pow5Factor(uint64_t value) {
184  const uint64_t m_inv_5 = 14757395258967641293u; // 5 * m_inv_5 = 1 (mod 2^64)
185  const uint64_t n_div_5 = 3689348814741910323u;  // #{ n | n = 0 (mod 2^64) } = 2^64 / 5
186  uint32_t count = 0;
187  for (;;) {
188    assert(value != 0);
189    value *= m_inv_5;
190    if (value > n_div_5)
191      break;
192    ++count;
193  }
194  return count;
195}
196
197// Returns true if value is divisible by 5^p.
198static inline bool multipleOfPowerOf5(const uint64_t value, const uint32_t p) {
199  // I tried a case distinction on p, but there was no performance difference.
200  return pow5Factor(value) >= p;
201}
202
203// Returns true if value is divisible by 2^p.
204static inline bool multipleOfPowerOf2(const uint64_t value, const uint32_t p) {
205  assert(value != 0);
206  assert(p < 64);
207  // __builtin_ctzll doesn't appear to be faster here.
208  return (value & ((1ull << p) - 1)) == 0;
209}
210
211// We need a 64x128-bit multiplication and a subsequent 128-bit shift.
212// Multiplication:
213//   The 64-bit factor is variable and passed in, the 128-bit factor comes
214//   from a lookup table. We know that the 64-bit factor only has 55
215//   significant bits (i.e., the 9 topmost bits are zeros). The 128-bit
216//   factor only has 124 significant bits (i.e., the 4 topmost bits are
217//   zeros).
218// Shift:
219//   In principle, the multiplication result requires 55 + 124 = 179 bits to
220//   represent. However, we then shift this value to the right by j, which is
221//   at least j >= 115, so the result is guaranteed to fit into 179 - 115 = 64
222//   bits. This means that we only need the topmost 64 significant bits of
223//   the 64x128-bit multiplication.
224//
225// There are several ways to do this:
226// 1. Best case: the compiler exposes a 128-bit type.
227//    We perform two 64x64-bit multiplications, add the higher 64 bits of the
228//    lower result to the higher result, and shift by j - 64 bits.
229//
230//    We explicitly cast from 64-bit to 128-bit, so the compiler can tell
231//    that these are only 64-bit inputs, and can map these to the best
232//    possible sequence of assembly instructions.
233//    x64 machines happen to have matching assembly instructions for
234//    64x64-bit multiplications and 128-bit shifts.
235//
236// 2. Second best case: the compiler exposes intrinsics for the x64 assembly
237//    instructions mentioned in 1.
238//
239// 3. We only have 64x64 bit instructions that return the lower 64 bits of
240//    the result, i.e., we have to use plain C.
241//    Our inputs are less than the full width, so we have three options:
242//    a. Ignore this fact and just implement the intrinsics manually.
243//    b. Split both into 31-bit pieces, which guarantees no internal overflow,
244//       but requires extra work upfront (unless we change the lookup table).
245//    c. Split only the first factor into 31-bit pieces, which also guarantees
246//       no internal overflow, but requires extra work since the intermediate
247//       results are not perfectly aligned.
248#if defined(HAS_UINT128)
249
250// Best case: use 128-bit type.
251static inline uint64_t mulShift64(const uint64_t m, const uint64_t* const mul, const int32_t j) {
252  const uint128_t b0 = ((uint128_t) m) * mul[0];
253  const uint128_t b2 = ((uint128_t) m) * mul[1];
254  return (uint64_t) (((b0 >> 64) + b2) >> (j - 64));
255}
256
257static inline uint64_t mulShiftAll64(const uint64_t m, const uint64_t* const mul, const int32_t j,
258  uint64_t* const vp, uint64_t* const vm, const uint32_t mmShift) {
259//  m <<= 2;
260//  uint128_t b0 = ((uint128_t) m) * mul[0]; // 0
261//  uint128_t b2 = ((uint128_t) m) * mul[1]; // 64
262//
263//  uint128_t hi = (b0 >> 64) + b2;
264//  uint128_t lo = b0 & 0xffffffffffffffffull;
265//  uint128_t factor = (((uint128_t) mul[1]) << 64) + mul[0];
266//  uint128_t vpLo = lo + (factor << 1);
267//  *vp = (uint64_t) ((hi + (vpLo >> 64)) >> (j - 64));
268//  uint128_t vmLo = lo - (factor << mmShift);
269//  *vm = (uint64_t) ((hi + (vmLo >> 64) - (((uint128_t) 1ull) << 64)) >> (j - 64));
270//  return (uint64_t) (hi >> (j - 64));
271  *vp = mulShift64(4 * m + 2, mul, j);
272  *vm = mulShift64(4 * m - 1 - mmShift, mul, j);
273  return mulShift64(4 * m, mul, j);
274}
275
276#elif defined(HAS_64_BIT_INTRINSICS)
277
278static inline uint64_t mulShift64(const uint64_t m, const uint64_t* const mul, const int32_t j) {
279  // m is maximum 55 bits
280  uint64_t high1;                                   // 128
281  const uint64_t low1 = umul128(m, mul[1], &high1); // 64
282  uint64_t high0;                                   // 64
283  umul128(m, mul[0], &high0);                       // 0
284  const uint64_t sum = high0 + low1;
285  if (sum < high0) {
286    ++high1; // overflow into high1
287  }
288  return shiftright128(sum, high1, j - 64);
289}
290
291static inline uint64_t mulShiftAll64(const uint64_t m, const uint64_t* const mul, const int32_t j,
292  uint64_t* const vp, uint64_t* const vm, const uint32_t mmShift) {
293  *vp = mulShift64(4 * m + 2, mul, j);
294  *vm = mulShift64(4 * m - 1 - mmShift, mul, j);
295  return mulShift64(4 * m, mul, j);
296}
297
298#else // !defined(HAS_UINT128) && !defined(HAS_64_BIT_INTRINSICS)
299
300static inline uint64_t mulShift64(const uint64_t m, const uint64_t* const mul, const int32_t j) {
301  // m is maximum 55 bits
302  uint64_t high1;                                   // 128
303  const uint64_t low1 = umul128(m, mul[1], &high1); // 64
304  uint64_t high0;                                   // 64
305  umul128(m, mul[0], &high0);                       // 0
306  const uint64_t sum = high0 + low1;
307  if (sum < high0) {
308    ++high1; // overflow into high1
309  }
310  return shiftright128(sum, high1, j - 64);
311}
312
313// This is faster if we don't have a 64x64->128-bit multiplication.
314static inline uint64_t mulShiftAll64(uint64_t m, const uint64_t* const mul, const int32_t j,
315  uint64_t* const vp, uint64_t* const vm, const uint32_t mmShift) {
316  m <<= 1;
317  // m is maximum 55 bits
318  uint64_t tmp;
319  const uint64_t lo = umul128(m, mul[0], &tmp);
320  uint64_t hi;
321  const uint64_t mid = tmp + umul128(m, mul[1], &hi);
322  hi += mid < tmp; // overflow into hi
323
324  const uint64_t lo2 = lo + mul[0];
325  const uint64_t mid2 = mid + mul[1] + (lo2 < lo);
326  const uint64_t hi2 = hi + (mid2 < mid);
327  *vp = shiftright128(mid2, hi2, (uint32_t) (j - 64 - 1));
328
329  if (mmShift == 1) {
330    const uint64_t lo3 = lo - mul[0];
331    const uint64_t mid3 = mid - mul[1] - (lo3 > lo);
332    const uint64_t hi3 = hi - (mid3 > mid);
333    *vm = shiftright128(mid3, hi3, (uint32_t) (j - 64 - 1));
334  } else {
335    const uint64_t lo3 = lo + lo;
336    const uint64_t mid3 = mid + mid + (lo3 < lo);
337    const uint64_t hi3 = hi + hi + (mid3 < mid);
338    const uint64_t lo4 = lo3 - mul[0];
339    const uint64_t mid4 = mid3 - mul[1] - (lo4 > lo3);
340    const uint64_t hi4 = hi3 - (mid4 > mid3);
341    *vm = shiftright128(mid4, hi4, (uint32_t) (j - 64));
342  }
343
344  return shiftright128(mid, hi, (uint32_t) (j - 64 - 1));
345}
346
347#endif // HAS_64_BIT_INTRINSICS
348
349#endif // RYU_D2S_INTRINSICS_H
350