rapidjson/internal/diyfp.h

0001 // Tencent is pleased to support the open source community by making RapidJSON available.
0002 //
0003 // Copyright (C) 2015 THL A29 Limited, a Tencent company, and Milo Yip.
0004 //
0005 // Licensed under the MIT License (the "License"); you may not use this file except
0006 // in compliance with the License. You may obtain a copy of the License at
0007 //
0008 // http://opensource.org/licenses/MIT
0009 //
0010 // Unless required by applicable law or agreed to in writing, software distributed
0011 // under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR
0012 // CONDITIONS OF ANY KIND, either express or implied. See the License for the
0013 // specific language governing permissions and limitations under the License.
0014
0015 // This is a C++ header-only implementation of Grisu2 algorithm from the publication:
0016 // Loitsch, Florian. "Printing floating-point numbers quickly and accurately with
0017 // integers." ACM Sigplan Notices 45.6 (2010): 233-243.
0018
0019 #ifndef RAPIDJSON_DIYFP_H_
0020 #define RAPIDJSON_DIYFP_H_
0021
0022 #include "../rapidjson.h"
0023 #include "clzll.h"
0024 #include <limits>
0025
0026 #if defined(_MSC_VER) && defined(_M_AMD64) && !defined(__INTEL_COMPILER)
0027 #include <intrin.h>
0028 #if !defined(_ARM64EC_)
0029 #pragma intrinsic(_umul128)
0030 #else
0031 #pragma comment(lib,"softintrin")
0032 #endif
0033 #endif
0034
0035 RAPIDJSON_NAMESPACE_BEGIN
0036 namespace internal {
0037
0038 #ifdef __GNUC__
0039 RAPIDJSON_DIAG_PUSH
0040 RAPIDJSON_DIAG_OFF(effc++)
0041 #endif
0042
0043 #ifdef __clang__
0044 RAPIDJSON_DIAG_PUSH
0045 RAPIDJSON_DIAG_OFF(padded)
0046 #endif
0047
0048 struct DiyFp {
0049     DiyFp() : f(), e() {}
0050
0051     DiyFp(uint64_t fp, int exp) : f(fp), e(exp) {}
0052
0053     explicit DiyFp(double d) {
0054         union {
0055             double d;
0056             uint64_t u64;
0057         } u = { d };
0058
0059         int biased_e = static_cast<int>((u.u64 & kDpExponentMask) >> kDpSignificandSize);
0060         uint64_t significand = (u.u64 & kDpSignificandMask);
0061         if (biased_e != 0) {
0062             f = significand + kDpHiddenBit;
0063             e = biased_e - kDpExponentBias;
0064         }
0065         else {
0066             f = significand;
0067             e = kDpMinExponent + 1;
0068         }
0069     }
0070
0071     DiyFp operator-(const DiyFp& rhs) const {
0072         return DiyFp(f - rhs.f, e);
0073     }
0074
0075     DiyFp operator*(const DiyFp& rhs) const {
0076 #if defined(_MSC_VER) && defined(_M_AMD64)
0077         uint64_t h;
0078         uint64_t l = _umul128(f, rhs.f, &h);
0079         if (l & (uint64_t(1) << 63)) // rounding
0080             h++;
0081         return DiyFp(h, e + rhs.e + 64);
0082 #elif defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)) && defined(__x86_64__)
0083         __extension__ typedef unsigned __int128 uint128;
0084         uint128 p = static_cast<uint128>(f) * static_cast<uint128>(rhs.f);
0085         uint64_t h = static_cast<uint64_t>(p >> 64);
0086         uint64_t l = static_cast<uint64_t>(p);
0087         if (l & (uint64_t(1) << 63)) // rounding
0088             h++;
0089         return DiyFp(h, e + rhs.e + 64);
0090 #else
0091         const uint64_t M32 = 0xFFFFFFFF;
0092         const uint64_t a = f >> 32;
0093         const uint64_t b = f & M32;
0094         const uint64_t c = rhs.f >> 32;
0095         const uint64_t d = rhs.f & M32;
0096         const uint64_t ac = a * c;
0097         const uint64_t bc = b * c;
0098         const uint64_t ad = a * d;
0099         const uint64_t bd = b * d;
0100         uint64_t tmp = (bd >> 32) + (ad & M32) + (bc & M32);
0101         tmp += 1U << 31;  /// mult_round
0102         return DiyFp(ac + (ad >> 32) + (bc >> 32) + (tmp >> 32), e + rhs.e + 64);
0103 #endif
0104     }
0105
0106     DiyFp Normalize() const {
0107         int s = static_cast<int>(clzll(f));
0108         return DiyFp(f << s, e - s);
0109     }
0110
0111     DiyFp NormalizeBoundary() const {
0112         DiyFp res = *this;
0113         while (!(res.f & (kDpHiddenBit << 1))) {
0114             res.f <<= 1;
0115             res.e--;
0116         }
0117         res.f <<= (kDiySignificandSize - kDpSignificandSize - 2);
0118         res.e = res.e - (kDiySignificandSize - kDpSignificandSize - 2);
0119         return res;
0120     }
0121
0122     void NormalizedBoundaries(DiyFp* minus, DiyFp* plus) const {
0123         DiyFp pl = DiyFp((f << 1) + 1, e - 1).NormalizeBoundary();
0124         DiyFp mi = (f == kDpHiddenBit) ? DiyFp((f << 2) - 1, e - 2) : DiyFp((f << 1) - 1, e - 1);
0125         mi.f <<= mi.e - pl.e;
0126         mi.e = pl.e;
0127         *plus = pl;
0128         *minus = mi;
0129     }
0130
0131     double ToDouble() const {
0132         union {
0133             double d;
0134             uint64_t u64;
0135         }u;
0136         RAPIDJSON_ASSERT(f <= kDpHiddenBit + kDpSignificandMask);
0137         if (e < kDpDenormalExponent) {
0138             // Underflow.
0139             return 0.0;
0140         }
0141         if (e >= kDpMaxExponent) {
0142             // Overflow.
0143             return std::numeric_limits<double>::infinity();
0144         }
0145         const uint64_t be = (e == kDpDenormalExponent && (f & kDpHiddenBit) == 0) ? 0 :
0146             static_cast<uint64_t>(e + kDpExponentBias);
0147         u.u64 = (f & kDpSignificandMask) | (be << kDpSignificandSize);
0148         return u.d;
0149     }
0150
0151     static const int kDiySignificandSize = 64;
0152     static const int kDpSignificandSize = 52;
0153     static const int kDpExponentBias = 0x3FF + kDpSignificandSize;
0154     static const int kDpMaxExponent = 0x7FF - kDpExponentBias;
0155     static const int kDpMinExponent = -kDpExponentBias;
0156     static const int kDpDenormalExponent = -kDpExponentBias + 1;
0157     static const uint64_t kDpExponentMask = RAPIDJSON_UINT64_C2(0x7FF00000, 0x00000000);
0158     static const uint64_t kDpSignificandMask = RAPIDJSON_UINT64_C2(0x000FFFFF, 0xFFFFFFFF);
0159     static const uint64_t kDpHiddenBit = RAPIDJSON_UINT64_C2(0x00100000, 0x00000000);
0160
0161     uint64_t f;
0162     int e;
0163 };
0164
0165 inline DiyFp GetCachedPowerByIndex(size_t index) {
0166     // 10^-348, 10^-340, ..., 10^340
0167     static const uint64_t kCachedPowers_F[] = {
0168         RAPIDJSON_UINT64_C2(0xfa8fd5a0, 0x081c0288), RAPIDJSON_UINT64_C2(0xbaaee17f, 0xa23ebf76),
0169         RAPIDJSON_UINT64_C2(0x8b16fb20, 0x3055ac76), RAPIDJSON_UINT64_C2(0xcf42894a, 0x5dce35ea),
0170         RAPIDJSON_UINT64_C2(0x9a6bb0aa, 0x55653b2d), RAPIDJSON_UINT64_C2(0xe61acf03, 0x3d1a45df),
0171         RAPIDJSON_UINT64_C2(0xab70fe17, 0xc79ac6ca), RAPIDJSON_UINT64_C2(0xff77b1fc, 0xbebcdc4f),
0172         RAPIDJSON_UINT64_C2(0xbe5691ef, 0x416bd60c), RAPIDJSON_UINT64_C2(0x8dd01fad, 0x907ffc3c),
0173         RAPIDJSON_UINT64_C2(0xd3515c28, 0x31559a83), RAPIDJSON_UINT64_C2(0x9d71ac8f, 0xada6c9b5),
0174         RAPIDJSON_UINT64_C2(0xea9c2277, 0x23ee8bcb), RAPIDJSON_UINT64_C2(0xaecc4991, 0x4078536d),
0175         RAPIDJSON_UINT64_C2(0x823c1279, 0x5db6ce57), RAPIDJSON_UINT64_C2(0xc2109436, 0x4dfb5637),
0176         RAPIDJSON_UINT64_C2(0x9096ea6f, 0x3848984f), RAPIDJSON_UINT64_C2(0xd77485cb, 0x25823ac7),
0177         RAPIDJSON_UINT64_C2(0xa086cfcd, 0x97bf97f4), RAPIDJSON_UINT64_C2(0xef340a98, 0x172aace5),
0178         RAPIDJSON_UINT64_C2(0xb23867fb, 0x2a35b28e), RAPIDJSON_UINT64_C2(0x84c8d4df, 0xd2c63f3b),
0179         RAPIDJSON_UINT64_C2(0xc5dd4427, 0x1ad3cdba), RAPIDJSON_UINT64_C2(0x936b9fce, 0xbb25c996),
0180         RAPIDJSON_UINT64_C2(0xdbac6c24, 0x7d62a584), RAPIDJSON_UINT64_C2(0xa3ab6658, 0x0d5fdaf6),
0181         RAPIDJSON_UINT64_C2(0xf3e2f893, 0xdec3f126), RAPIDJSON_UINT64_C2(0xb5b5ada8, 0xaaff80b8),
0182         RAPIDJSON_UINT64_C2(0x87625f05, 0x6c7c4a8b), RAPIDJSON_UINT64_C2(0xc9bcff60, 0x34c13053),
0183         RAPIDJSON_UINT64_C2(0x964e858c, 0x91ba2655), RAPIDJSON_UINT64_C2(0xdff97724, 0x70297ebd),
0184         RAPIDJSON_UINT64_C2(0xa6dfbd9f, 0xb8e5b88f), RAPIDJSON_UINT64_C2(0xf8a95fcf, 0x88747d94),
0185         RAPIDJSON_UINT64_C2(0xb9447093, 0x8fa89bcf), RAPIDJSON_UINT64_C2(0x8a08f0f8, 0xbf0f156b),
0186         RAPIDJSON_UINT64_C2(0xcdb02555, 0x653131b6), RAPIDJSON_UINT64_C2(0x993fe2c6, 0xd07b7fac),
0187         RAPIDJSON_UINT64_C2(0xe45c10c4, 0x2a2b3b06), RAPIDJSON_UINT64_C2(0xaa242499, 0x697392d3),
0188         RAPIDJSON_UINT64_C2(0xfd87b5f2, 0x8300ca0e), RAPIDJSON_UINT64_C2(0xbce50864, 0x92111aeb),
0189         RAPIDJSON_UINT64_C2(0x8cbccc09, 0x6f5088cc), RAPIDJSON_UINT64_C2(0xd1b71758, 0xe219652c),
0190         RAPIDJSON_UINT64_C2(0x9c400000, 0x00000000), RAPIDJSON_UINT64_C2(0xe8d4a510, 0x00000000),
0191         RAPIDJSON_UINT64_C2(0xad78ebc5, 0xac620000), RAPIDJSON_UINT64_C2(0x813f3978, 0xf8940984),
0192         RAPIDJSON_UINT64_C2(0xc097ce7b, 0xc90715b3), RAPIDJSON_UINT64_C2(0x8f7e32ce, 0x7bea5c70),
0193         RAPIDJSON_UINT64_C2(0xd5d238a4, 0xabe98068), RAPIDJSON_UINT64_C2(0x9f4f2726, 0x179a2245),
0194         RAPIDJSON_UINT64_C2(0xed63a231, 0xd4c4fb27), RAPIDJSON_UINT64_C2(0xb0de6538, 0x8cc8ada8),
0195         RAPIDJSON_UINT64_C2(0x83c7088e, 0x1aab65db), RAPIDJSON_UINT64_C2(0xc45d1df9, 0x42711d9a),
0196         RAPIDJSON_UINT64_C2(0x924d692c, 0xa61be758), RAPIDJSON_UINT64_C2(0xda01ee64, 0x1a708dea),
0197         RAPIDJSON_UINT64_C2(0xa26da399, 0x9aef774a), RAPIDJSON_UINT64_C2(0xf209787b, 0xb47d6b85),
0198         RAPIDJSON_UINT64_C2(0xb454e4a1, 0x79dd1877), RAPIDJSON_UINT64_C2(0x865b8692, 0x5b9bc5c2),
0199         RAPIDJSON_UINT64_C2(0xc83553c5, 0xc8965d3d), RAPIDJSON_UINT64_C2(0x952ab45c, 0xfa97a0b3),
0200         RAPIDJSON_UINT64_C2(0xde469fbd, 0x99a05fe3), RAPIDJSON_UINT64_C2(0xa59bc234, 0xdb398c25),
0201         RAPIDJSON_UINT64_C2(0xf6c69a72, 0xa3989f5c), RAPIDJSON_UINT64_C2(0xb7dcbf53, 0x54e9bece),
0202         RAPIDJSON_UINT64_C2(0x88fcf317, 0xf22241e2), RAPIDJSON_UINT64_C2(0xcc20ce9b, 0xd35c78a5),
0203         RAPIDJSON_UINT64_C2(0x98165af3, 0x7b2153df), RAPIDJSON_UINT64_C2(0xe2a0b5dc, 0x971f303a),
0204         RAPIDJSON_UINT64_C2(0xa8d9d153, 0x5ce3b396), RAPIDJSON_UINT64_C2(0xfb9b7cd9, 0xa4a7443c),
0205         RAPIDJSON_UINT64_C2(0xbb764c4c, 0xa7a44410), RAPIDJSON_UINT64_C2(0x8bab8eef, 0xb6409c1a),
0206         RAPIDJSON_UINT64_C2(0xd01fef10, 0xa657842c), RAPIDJSON_UINT64_C2(0x9b10a4e5, 0xe9913129),
0207         RAPIDJSON_UINT64_C2(0xe7109bfb, 0xa19c0c9d), RAPIDJSON_UINT64_C2(0xac2820d9, 0x623bf429),
0208         RAPIDJSON_UINT64_C2(0x80444b5e, 0x7aa7cf85), RAPIDJSON_UINT64_C2(0xbf21e440, 0x03acdd2d),
0209         RAPIDJSON_UINT64_C2(0x8e679c2f, 0x5e44ff8f), RAPIDJSON_UINT64_C2(0xd433179d, 0x9c8cb841),
0210         RAPIDJSON_UINT64_C2(0x9e19db92, 0xb4e31ba9), RAPIDJSON_UINT64_C2(0xeb96bf6e, 0xbadf77d9),
0211         RAPIDJSON_UINT64_C2(0xaf87023b, 0x9bf0ee6b)
0212     };
0213     static const int16_t kCachedPowers_E[] = {
0214         -1220, -1193, -1166, -1140, -1113, -1087, -1060, -1034, -1007,  -980,
0215         -954,  -927,  -901,  -874,  -847,  -821,  -794,  -768,  -741,  -715,
0216         -688,  -661,  -635,  -608,  -582,  -555,  -529,  -502,  -475,  -449,
0217         -422,  -396,  -369,  -343,  -316,  -289,  -263,  -236,  -210,  -183,
0218         -157,  -130,  -103,   -77,   -50,   -24,     3,    30,    56,    83,
0219         109,   136,   162,   189,   216,   242,   269,   295,   322,   348,
0220         375,   402,   428,   455,   481,   508,   534,   561,   588,   614,
0221         641,   667,   694,   720,   747,   774,   800,   827,   853,   880,
0222         907,   933,   960,   986,  1013,  1039,  1066
0223     };
0224     RAPIDJSON_ASSERT(index < 87);
0225     return DiyFp(kCachedPowers_F[index], kCachedPowers_E[index]);
0226 }
0227
0228 inline DiyFp GetCachedPower(int e, int* K) {
0229
0230     //int k = static_cast<int>(ceil((-61 - e) * 0.30102999566398114)) + 374;
0231     double dk = (-61 - e) * 0.30102999566398114 + 347;  // dk must be positive, so can do ceiling in positive
0232     int k = static_cast<int>(dk);
0233     if (dk - k > 0.0)
0234         k++;
0235
0236     unsigned index = static_cast<unsigned>((k >> 3) + 1);
0237     *K = -(-348 + static_cast<int>(index << 3));    // decimal exponent no need lookup table
0238
0239     return GetCachedPowerByIndex(index);
0240 }
0241
0242 inline DiyFp GetCachedPower10(int exp, int *outExp) {
0243     RAPIDJSON_ASSERT(exp >= -348);
0244     unsigned index = static_cast<unsigned>(exp + 348) / 8u;
0245     *outExp = -348 + static_cast<int>(index) * 8;
0246     return GetCachedPowerByIndex(index);
0247 }
0248
0249 #ifdef __GNUC__
0250 RAPIDJSON_DIAG_POP
0251 #endif
0252
0253 #ifdef __clang__
0254 RAPIDJSON_DIAG_POP
0255 RAPIDJSON_DIAG_OFF(padded)
0256 #endif
0257
0258 } // namespace internal
0259 RAPIDJSON_NAMESPACE_END
0260
0261 #endif // RAPIDJSON_DIYFP_H_