diyfp.h
1 // Tencent is pleased to support the open source community by making RapidJSON available.
2 //
3 // Copyright (C) 2015 THL A29 Limited, a Tencent company, and Milo Yip. All rights reserved.
4 //
5 // Licensed under the MIT License (the "License"); you may not use this file except
6 // in compliance with the License. You may obtain a copy of the License at
7 //
8 // http://opensource.org/licenses/MIT
9 //
10 // Unless required by applicable law or agreed to in writing, software distributed
11 // under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR
12 // CONDITIONS OF ANY KIND, either express or implied. See the License for the
13 // specific language governing permissions and limitations under the License.
14 
15 // This is a C++ header-only implementation of Grisu2 algorithm from the publication:
16 // Loitsch, Florian. "Printing floating-point numbers quickly and accurately with
17 // integers." ACM Sigplan Notices 45.6 (2010): 233-243.
18 
19 #ifndef RAPIDJSON_DIYFP_H_
20 #define RAPIDJSON_DIYFP_H_
21 
22 #include "../rapidjson.h"
23 
24 #if defined(_MSC_VER) && defined(_M_AMD64)
25 #include <intrin.h>
26 #pragma intrinsic(_BitScanReverse64)
27 #endif
28 
29 RAPIDJSON_NAMESPACE_BEGIN
30 namespace internal {
31 
32 #ifdef __GNUC__
33 RAPIDJSON_DIAG_PUSH
34 RAPIDJSON_DIAG_OFF(effc++)
35 #endif
36 
37 struct DiyFp {
38  DiyFp() {}
39 
40  DiyFp(uint64_t fp, int exp) : f(fp), e(exp) {}
41 
42  explicit DiyFp(double d) {
43  union {
44  double d;
45  uint64_t u64;
46  } u = { d };
47 
48  int biased_e = static_cast<int>((u.u64 & kDpExponentMask) >> kDpSignificandSize);
49  uint64_t significand = (u.u64 & kDpSignificandMask);
50  if (biased_e != 0) {
51  f = significand + kDpHiddenBit;
52  e = biased_e - kDpExponentBias;
53  }
54  else {
55  f = significand;
56  e = kDpMinExponent + 1;
57  }
58  }
59 
60  DiyFp operator-(const DiyFp& rhs) const {
61  return DiyFp(f - rhs.f, e);
62  }
63 
64  DiyFp operator*(const DiyFp& rhs) const {
65 #if defined(_MSC_VER) && defined(_M_AMD64)
66  uint64_t h;
67  uint64_t l = _umul128(f, rhs.f, &h);
68  if (l & (uint64_t(1) << 63)) // rounding
69  h++;
70  return DiyFp(h, e + rhs.e + 64);
71 #elif (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)) && defined(__x86_64__)
72  __extension__ typedef unsigned __int128 uint128;
73  uint128 p = static_cast<uint128>(f) * static_cast<uint128>(rhs.f);
74  uint64_t h = static_cast<uint64_t>(p >> 64);
75  uint64_t l = static_cast<uint64_t>(p);
76  if (l & (uint64_t(1) << 63)) // rounding
77  h++;
78  return DiyFp(h, e + rhs.e + 64);
79 #else
80  const uint64_t M32 = 0xFFFFFFFF;
81  const uint64_t a = f >> 32;
82  const uint64_t b = f & M32;
83  const uint64_t c = rhs.f >> 32;
84  const uint64_t d = rhs.f & M32;
85  const uint64_t ac = a * c;
86  const uint64_t bc = b * c;
87  const uint64_t ad = a * d;
88  const uint64_t bd = b * d;
89  uint64_t tmp = (bd >> 32) + (ad & M32) + (bc & M32);
90  tmp += 1U << 31; /// mult_round
91  return DiyFp(ac + (ad >> 32) + (bc >> 32) + (tmp >> 32), e + rhs.e + 64);
92 #endif
93  }
94 
95  DiyFp Normalize() const {
96 #if defined(_MSC_VER) && defined(_M_AMD64)
97  unsigned long index;
98  _BitScanReverse64(&index, f);
99  return DiyFp(f << (63 - index), e - (63 - index));
100 #elif defined(__GNUC__) && __GNUC__ >= 4
101  int s = __builtin_clzll(f);
102  return DiyFp(f << s, e - s);
103 #else
104  DiyFp res = *this;
105  while (!(res.f & (static_cast<uint64_t>(1) << 63))) {
106  res.f <<= 1;
107  res.e--;
108  }
109  return res;
110 #endif
111  }
112 
113  DiyFp NormalizeBoundary() const {
114  DiyFp res = *this;
115  while (!(res.f & (kDpHiddenBit << 1))) {
116  res.f <<= 1;
117  res.e--;
118  }
119  res.f <<= (kDiySignificandSize - kDpSignificandSize - 2);
120  res.e = res.e - (kDiySignificandSize - kDpSignificandSize - 2);
121  return res;
122  }
123 
124  void NormalizedBoundaries(DiyFp* minus, DiyFp* plus) const {
125  DiyFp pl = DiyFp((f << 1) + 1, e - 1).NormalizeBoundary();
126  DiyFp mi = (f == kDpHiddenBit) ? DiyFp((f << 2) - 1, e - 2) : DiyFp((f << 1) - 1, e - 1);
127  mi.f <<= mi.e - pl.e;
128  mi.e = pl.e;
129  *plus = pl;
130  *minus = mi;
131  }
132 
133  double ToDouble() const {
134  union {
135  double d;
136  uint64_t u64;
137  }u;
138  const uint64_t be = (e == kDpDenormalExponent && (f & kDpHiddenBit) == 0) ? 0 :
139  static_cast<uint64_t>(e + kDpExponentBias);
140  u.u64 = (f & kDpSignificandMask) | (be << kDpSignificandSize);
141  return u.d;
142  }
143 
144  static const int kDiySignificandSize = 64;
145  static const int kDpSignificandSize = 52;
146  static const int kDpExponentBias = 0x3FF + kDpSignificandSize;
147  static const int kDpMaxExponent = 0x7FF - kDpExponentBias;
148  static const int kDpMinExponent = -kDpExponentBias;
149  static const int kDpDenormalExponent = -kDpExponentBias + 1;
150  static const uint64_t kDpExponentMask = RAPIDJSON_UINT64_C2(0x7FF00000, 0x00000000);
151  static const uint64_t kDpSignificandMask = RAPIDJSON_UINT64_C2(0x000FFFFF, 0xFFFFFFFF);
152  static const uint64_t kDpHiddenBit = RAPIDJSON_UINT64_C2(0x00100000, 0x00000000);
153 
154  uint64_t f;
155  int e;
156 };
157 
158 inline DiyFp GetCachedPowerByIndex(size_t index) {
159  // 10^-348, 10^-340, ..., 10^340
160  static const uint64_t kCachedPowers_F[] = {
161  RAPIDJSON_UINT64_C2(0xfa8fd5a0, 0x081c0288), RAPIDJSON_UINT64_C2(0xbaaee17f, 0xa23ebf76),
162  RAPIDJSON_UINT64_C2(0x8b16fb20, 0x3055ac76), RAPIDJSON_UINT64_C2(0xcf42894a, 0x5dce35ea),
163  RAPIDJSON_UINT64_C2(0x9a6bb0aa, 0x55653b2d), RAPIDJSON_UINT64_C2(0xe61acf03, 0x3d1a45df),
164  RAPIDJSON_UINT64_C2(0xab70fe17, 0xc79ac6ca), RAPIDJSON_UINT64_C2(0xff77b1fc, 0xbebcdc4f),
165  RAPIDJSON_UINT64_C2(0xbe5691ef, 0x416bd60c), RAPIDJSON_UINT64_C2(0x8dd01fad, 0x907ffc3c),
166  RAPIDJSON_UINT64_C2(0xd3515c28, 0x31559a83), RAPIDJSON_UINT64_C2(0x9d71ac8f, 0xada6c9b5),
167  RAPIDJSON_UINT64_C2(0xea9c2277, 0x23ee8bcb), RAPIDJSON_UINT64_C2(0xaecc4991, 0x4078536d),
168  RAPIDJSON_UINT64_C2(0x823c1279, 0x5db6ce57), RAPIDJSON_UINT64_C2(0xc2109436, 0x4dfb5637),
169  RAPIDJSON_UINT64_C2(0x9096ea6f, 0x3848984f), RAPIDJSON_UINT64_C2(0xd77485cb, 0x25823ac7),
170  RAPIDJSON_UINT64_C2(0xa086cfcd, 0x97bf97f4), RAPIDJSON_UINT64_C2(0xef340a98, 0x172aace5),
171  RAPIDJSON_UINT64_C2(0xb23867fb, 0x2a35b28e), RAPIDJSON_UINT64_C2(0x84c8d4df, 0xd2c63f3b),
172  RAPIDJSON_UINT64_C2(0xc5dd4427, 0x1ad3cdba), RAPIDJSON_UINT64_C2(0x936b9fce, 0xbb25c996),
173  RAPIDJSON_UINT64_C2(0xdbac6c24, 0x7d62a584), RAPIDJSON_UINT64_C2(0xa3ab6658, 0x0d5fdaf6),
174  RAPIDJSON_UINT64_C2(0xf3e2f893, 0xdec3f126), RAPIDJSON_UINT64_C2(0xb5b5ada8, 0xaaff80b8),
175  RAPIDJSON_UINT64_C2(0x87625f05, 0x6c7c4a8b), RAPIDJSON_UINT64_C2(0xc9bcff60, 0x34c13053),
176  RAPIDJSON_UINT64_C2(0x964e858c, 0x91ba2655), RAPIDJSON_UINT64_C2(0xdff97724, 0x70297ebd),
177  RAPIDJSON_UINT64_C2(0xa6dfbd9f, 0xb8e5b88f), RAPIDJSON_UINT64_C2(0xf8a95fcf, 0x88747d94),
178  RAPIDJSON_UINT64_C2(0xb9447093, 0x8fa89bcf), RAPIDJSON_UINT64_C2(0x8a08f0f8, 0xbf0f156b),
179  RAPIDJSON_UINT64_C2(0xcdb02555, 0x653131b6), RAPIDJSON_UINT64_C2(0x993fe2c6, 0xd07b7fac),
180  RAPIDJSON_UINT64_C2(0xe45c10c4, 0x2a2b3b06), RAPIDJSON_UINT64_C2(0xaa242499, 0x697392d3),
181  RAPIDJSON_UINT64_C2(0xfd87b5f2, 0x8300ca0e), RAPIDJSON_UINT64_C2(0xbce50864, 0x92111aeb),
182  RAPIDJSON_UINT64_C2(0x8cbccc09, 0x6f5088cc), RAPIDJSON_UINT64_C2(0xd1b71758, 0xe219652c),
183  RAPIDJSON_UINT64_C2(0x9c400000, 0x00000000), RAPIDJSON_UINT64_C2(0xe8d4a510, 0x00000000),
184  RAPIDJSON_UINT64_C2(0xad78ebc5, 0xac620000), RAPIDJSON_UINT64_C2(0x813f3978, 0xf8940984),
185  RAPIDJSON_UINT64_C2(0xc097ce7b, 0xc90715b3), RAPIDJSON_UINT64_C2(0x8f7e32ce, 0x7bea5c70),
186  RAPIDJSON_UINT64_C2(0xd5d238a4, 0xabe98068), RAPIDJSON_UINT64_C2(0x9f4f2726, 0x179a2245),
187  RAPIDJSON_UINT64_C2(0xed63a231, 0xd4c4fb27), RAPIDJSON_UINT64_C2(0xb0de6538, 0x8cc8ada8),
188  RAPIDJSON_UINT64_C2(0x83c7088e, 0x1aab65db), RAPIDJSON_UINT64_C2(0xc45d1df9, 0x42711d9a),
189  RAPIDJSON_UINT64_C2(0x924d692c, 0xa61be758), RAPIDJSON_UINT64_C2(0xda01ee64, 0x1a708dea),
190  RAPIDJSON_UINT64_C2(0xa26da399, 0x9aef774a), RAPIDJSON_UINT64_C2(0xf209787b, 0xb47d6b85),
191  RAPIDJSON_UINT64_C2(0xb454e4a1, 0x79dd1877), RAPIDJSON_UINT64_C2(0x865b8692, 0x5b9bc5c2),
192  RAPIDJSON_UINT64_C2(0xc83553c5, 0xc8965d3d), RAPIDJSON_UINT64_C2(0x952ab45c, 0xfa97a0b3),
193  RAPIDJSON_UINT64_C2(0xde469fbd, 0x99a05fe3), RAPIDJSON_UINT64_C2(0xa59bc234, 0xdb398c25),
194  RAPIDJSON_UINT64_C2(0xf6c69a72, 0xa3989f5c), RAPIDJSON_UINT64_C2(0xb7dcbf53, 0x54e9bece),
195  RAPIDJSON_UINT64_C2(0x88fcf317, 0xf22241e2), RAPIDJSON_UINT64_C2(0xcc20ce9b, 0xd35c78a5),
196  RAPIDJSON_UINT64_C2(0x98165af3, 0x7b2153df), RAPIDJSON_UINT64_C2(0xe2a0b5dc, 0x971f303a),
197  RAPIDJSON_UINT64_C2(0xa8d9d153, 0x5ce3b396), RAPIDJSON_UINT64_C2(0xfb9b7cd9, 0xa4a7443c),
198  RAPIDJSON_UINT64_C2(0xbb764c4c, 0xa7a44410), RAPIDJSON_UINT64_C2(0x8bab8eef, 0xb6409c1a),
199  RAPIDJSON_UINT64_C2(0xd01fef10, 0xa657842c), RAPIDJSON_UINT64_C2(0x9b10a4e5, 0xe9913129),
200  RAPIDJSON_UINT64_C2(0xe7109bfb, 0xa19c0c9d), RAPIDJSON_UINT64_C2(0xac2820d9, 0x623bf429),
201  RAPIDJSON_UINT64_C2(0x80444b5e, 0x7aa7cf85), RAPIDJSON_UINT64_C2(0xbf21e440, 0x03acdd2d),
202  RAPIDJSON_UINT64_C2(0x8e679c2f, 0x5e44ff8f), RAPIDJSON_UINT64_C2(0xd433179d, 0x9c8cb841),
203  RAPIDJSON_UINT64_C2(0x9e19db92, 0xb4e31ba9), RAPIDJSON_UINT64_C2(0xeb96bf6e, 0xbadf77d9),
204  RAPIDJSON_UINT64_C2(0xaf87023b, 0x9bf0ee6b)
205  };
206  static const int16_t kCachedPowers_E[] = {
207  -1220, -1193, -1166, -1140, -1113, -1087, -1060, -1034, -1007, -980,
208  -954, -927, -901, -874, -847, -821, -794, -768, -741, -715,
209  -688, -661, -635, -608, -582, -555, -529, -502, -475, -449,
210  -422, -396, -369, -343, -316, -289, -263, -236, -210, -183,
211  -157, -130, -103, -77, -50, -24, 3, 30, 56, 83,
212  109, 136, 162, 189, 216, 242, 269, 295, 322, 348,
213  375, 402, 428, 455, 481, 508, 534, 561, 588, 614,
214  641, 667, 694, 720, 747, 774, 800, 827, 853, 880,
215  907, 933, 960, 986, 1013, 1039, 1066
216  };
217  return DiyFp(kCachedPowers_F[index], kCachedPowers_E[index]);
218 }
219 
220 inline DiyFp GetCachedPower(int e, int* K) {
221 
222  //int k = static_cast<int>(ceil((-61 - e) * 0.30102999566398114)) + 374;
223  double dk = (-61 - e) * 0.30102999566398114 + 347; // dk must be positive, so can do ceiling in positive
224  int k = static_cast<int>(dk);
225  if (dk - k > 0.0)
226  k++;
227 
228  unsigned index = static_cast<unsigned>((k >> 3) + 1);
229  *K = -(-348 + static_cast<int>(index << 3)); // decimal exponent no need lookup table
230 
231  return GetCachedPowerByIndex(index);
232 }
233 
234 inline DiyFp GetCachedPower10(int exp, int *outExp) {
235  unsigned index = (exp + 348) / 8;
236  *outExp = -348 + index * 8;
237  return GetCachedPowerByIndex(index);
238  }
239 
240 #ifdef __GNUC__
241 RAPIDJSON_DIAG_POP
242 #endif
243 
244 } // namespace internal
245 RAPIDJSON_NAMESPACE_END
246 
247 #endif // RAPIDJSON_DIYFP_H_
#define RAPIDJSON_UINT64_C2(high32, low32)
Construct a 64-bit literal by a pair of 32-bit integer.
Definition: rapidjson.h:261