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