diyfp.h

00001 // Tencent is pleased to support the open source community by making RapidJSON available.
00002 // 
00003 // Copyright (C) 2015 THL A29 Limited, a Tencent company, and Milo Yip. All rights reserved.
00004 //
00005 // Licensed under the MIT License (the "License"); you may not use this file except
00006 // in compliance with the License. You may obtain a copy of the License at
00007 //
00008 // http://opensource.org/licenses/MIT
00009 //
00010 // Unless required by applicable law or agreed to in writing, software distributed 
00011 // under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR 
00012 // CONDITIONS OF ANY KIND, either express or implied. See the License for the 
00013 // specific language governing permissions and limitations under the License.
00014 
00015 // This is a C++ header-only implementation of Grisu2 algorithm from the publication:
00016 // Loitsch, Florian. "Printing floating-point numbers quickly and accurately with
00017 // integers." ACM Sigplan Notices 45.6 (2010): 233-243.
00018 
00019 #ifndef RAPIDJSON_DIYFP_H_
00020 #define RAPIDJSON_DIYFP_H_
00021 
00022 #include "../rapidjson.h"
00023 
00024 #if defined(_MSC_VER) && defined(_M_AMD64)
00025 #include <intrin.h>
00026 #pragma intrinsic(_BitScanReverse64)
00027 #pragma intrinsic(_umul128)
00028 #endif
00029 
00030 RAPIDJSON_NAMESPACE_BEGIN
00031 namespace internal {
00032 
00033 #ifdef __GNUC__
00034 RAPIDJSON_DIAG_PUSH
00035 RAPIDJSON_DIAG_OFF(effc++)
00036 #endif
00037 
00038 #ifdef __clang__
00039 RAPIDJSON_DIAG_PUSH
00040 RAPIDJSON_DIAG_OFF(padded)
00041 #endif
00042 
00043 struct DiyFp {
00044     DiyFp() : f(), e() {}
00045 
00046     DiyFp(uint64_t fp, int exp) : f(fp), e(exp) {}
00047 
00048     explicit DiyFp(double d) {
00049         union {
00050             double d;
00051             uint64_t u64;
00052         } u = { d };
00053 
00054         int biased_e = static_cast<int>((u.u64 & kDpExponentMask) >> kDpSignificandSize);
00055         uint64_t significand = (u.u64 & kDpSignificandMask);
00056         if (biased_e != 0) {
00057             f = significand + kDpHiddenBit;
00058             e = biased_e - kDpExponentBias;
00059         } 
00060         else {
00061             f = significand;
00062             e = kDpMinExponent + 1;
00063         }
00064     }
00065 
00066     DiyFp operator-(const DiyFp& rhs) const {
00067         return DiyFp(f - rhs.f, e);
00068     }
00069 
00070     DiyFp operator*(const DiyFp& rhs) const {
00071 #if defined(_MSC_VER) && defined(_M_AMD64)
00072         uint64_t h;
00073         uint64_t l = _umul128(f, rhs.f, &h);
00074         if (l & (uint64_t(1) << 63)) // rounding
00075             h++;
00076         return DiyFp(h, e + rhs.e + 64);
00077 #elif (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)) && defined(__x86_64__)
00078         __extension__ typedef unsigned __int128 uint128;
00079         uint128 p = static_cast<uint128>(f) * static_cast<uint128>(rhs.f);
00080         uint64_t h = static_cast<uint64_t>(p >> 64);
00081         uint64_t l = static_cast<uint64_t>(p);
00082         if (l & (uint64_t(1) << 63)) // rounding
00083             h++;
00084         return DiyFp(h, e + rhs.e + 64);
00085 #else
00086         const uint64_t M32 = 0xFFFFFFFF;
00087         const uint64_t a = f >> 32;
00088         const uint64_t b = f & M32;
00089         const uint64_t c = rhs.f >> 32;
00090         const uint64_t d = rhs.f & M32;
00091         const uint64_t ac = a * c;
00092         const uint64_t bc = b * c;
00093         const uint64_t ad = a * d;
00094         const uint64_t bd = b * d;
00095         uint64_t tmp = (bd >> 32) + (ad & M32) + (bc & M32);
00096         tmp += 1U << 31;  /// mult_round
00097         return DiyFp(ac + (ad >> 32) + (bc >> 32) + (tmp >> 32), e + rhs.e + 64);
00098 #endif
00099     }
00100 
00101     DiyFp Normalize() const {
00102 #if defined(_MSC_VER) && defined(_M_AMD64)
00103         unsigned long index;
00104         _BitScanReverse64(&index, f);
00105         return DiyFp(f << (63 - index), e - (63 - index));
00106 #elif defined(__GNUC__) && __GNUC__ >= 4
00107         int s = __builtin_clzll(f);
00108         return DiyFp(f << s, e - s);
00109 #else
00110         DiyFp res = *this;
00111         while (!(res.f & (static_cast<uint64_t>(1) << 63))) {
00112             res.f <<= 1;
00113             res.e--;
00114         }
00115         return res;
00116 #endif
00117     }
00118 
00119     DiyFp NormalizeBoundary() const {
00120         DiyFp res = *this;
00121         while (!(res.f & (kDpHiddenBit << 1))) {
00122             res.f <<= 1;
00123             res.e--;
00124         }
00125         res.f <<= (kDiySignificandSize - kDpSignificandSize - 2);
00126         res.e = res.e - (kDiySignificandSize - kDpSignificandSize - 2);
00127         return res;
00128     }
00129 
00130     void NormalizedBoundaries(DiyFp* minus, DiyFp* plus) const {
00131         DiyFp pl = DiyFp((f << 1) + 1, e - 1).NormalizeBoundary();
00132         DiyFp mi = (f == kDpHiddenBit) ? DiyFp((f << 2) - 1, e - 2) : DiyFp((f << 1) - 1, e - 1);
00133         mi.f <<= mi.e - pl.e;
00134         mi.e = pl.e;
00135         *plus = pl;
00136         *minus = mi;
00137     }
00138 
00139     double ToDouble() const {
00140         union {
00141             double d;
00142             uint64_t u64;
00143         }u;
00144         const uint64_t be = (e == kDpDenormalExponent && (f & kDpHiddenBit) == 0) ? 0 : 
00145             static_cast<uint64_t>(e + kDpExponentBias);
00146         u.u64 = (f & kDpSignificandMask) | (be << kDpSignificandSize);
00147         return u.d;
00148     }
00149 
00150     static const int kDiySignificandSize = 64;
00151     static const int kDpSignificandSize = 52;
00152     static const int kDpExponentBias = 0x3FF + kDpSignificandSize;
00153     static const int kDpMaxExponent = 0x7FF - kDpExponentBias;
00154     static const int kDpMinExponent = -kDpExponentBias;
00155     static const int kDpDenormalExponent = -kDpExponentBias + 1;
00156     static const uint64_t kDpExponentMask = RAPIDJSON_UINT64_C2(0x7FF00000, 0x00000000);
00157     static const uint64_t kDpSignificandMask = RAPIDJSON_UINT64_C2(0x000FFFFF, 0xFFFFFFFF);
00158     static const uint64_t kDpHiddenBit = RAPIDJSON_UINT64_C2(0x00100000, 0x00000000);
00159 
00160     uint64_t f;
00161     int e;
00162 };
00163 
00164 inline DiyFp GetCachedPowerByIndex(size_t index) {
00165     // 10^-348, 10^-340, ..., 10^340
00166     static const uint64_t kCachedPowers_F[] = {
00167         RAPIDJSON_UINT64_C2(0xfa8fd5a0, 0x081c0288), RAPIDJSON_UINT64_C2(0xbaaee17f, 0xa23ebf76),
00168         RAPIDJSON_UINT64_C2(0x8b16fb20, 0x3055ac76), RAPIDJSON_UINT64_C2(0xcf42894a, 0x5dce35ea),
00169         RAPIDJSON_UINT64_C2(0x9a6bb0aa, 0x55653b2d), RAPIDJSON_UINT64_C2(0xe61acf03, 0x3d1a45df),
00170         RAPIDJSON_UINT64_C2(0xab70fe17, 0xc79ac6ca), RAPIDJSON_UINT64_C2(0xff77b1fc, 0xbebcdc4f),
00171         RAPIDJSON_UINT64_C2(0xbe5691ef, 0x416bd60c), RAPIDJSON_UINT64_C2(0x8dd01fad, 0x907ffc3c),
00172         RAPIDJSON_UINT64_C2(0xd3515c28, 0x31559a83), RAPIDJSON_UINT64_C2(0x9d71ac8f, 0xada6c9b5),
00173         RAPIDJSON_UINT64_C2(0xea9c2277, 0x23ee8bcb), RAPIDJSON_UINT64_C2(0xaecc4991, 0x4078536d),
00174         RAPIDJSON_UINT64_C2(0x823c1279, 0x5db6ce57), RAPIDJSON_UINT64_C2(0xc2109436, 0x4dfb5637),
00175         RAPIDJSON_UINT64_C2(0x9096ea6f, 0x3848984f), RAPIDJSON_UINT64_C2(0xd77485cb, 0x25823ac7),
00176         RAPIDJSON_UINT64_C2(0xa086cfcd, 0x97bf97f4), RAPIDJSON_UINT64_C2(0xef340a98, 0x172aace5),
00177         RAPIDJSON_UINT64_C2(0xb23867fb, 0x2a35b28e), RAPIDJSON_UINT64_C2(0x84c8d4df, 0xd2c63f3b),
00178         RAPIDJSON_UINT64_C2(0xc5dd4427, 0x1ad3cdba), RAPIDJSON_UINT64_C2(0x936b9fce, 0xbb25c996),
00179         RAPIDJSON_UINT64_C2(0xdbac6c24, 0x7d62a584), RAPIDJSON_UINT64_C2(0xa3ab6658, 0x0d5fdaf6),
00180         RAPIDJSON_UINT64_C2(0xf3e2f893, 0xdec3f126), RAPIDJSON_UINT64_C2(0xb5b5ada8, 0xaaff80b8),
00181         RAPIDJSON_UINT64_C2(0x87625f05, 0x6c7c4a8b), RAPIDJSON_UINT64_C2(0xc9bcff60, 0x34c13053),
00182         RAPIDJSON_UINT64_C2(0x964e858c, 0x91ba2655), RAPIDJSON_UINT64_C2(0xdff97724, 0x70297ebd),
00183         RAPIDJSON_UINT64_C2(0xa6dfbd9f, 0xb8e5b88f), RAPIDJSON_UINT64_C2(0xf8a95fcf, 0x88747d94),
00184         RAPIDJSON_UINT64_C2(0xb9447093, 0x8fa89bcf), RAPIDJSON_UINT64_C2(0x8a08f0f8, 0xbf0f156b),
00185         RAPIDJSON_UINT64_C2(0xcdb02555, 0x653131b6), RAPIDJSON_UINT64_C2(0x993fe2c6, 0xd07b7fac),
00186         RAPIDJSON_UINT64_C2(0xe45c10c4, 0x2a2b3b06), RAPIDJSON_UINT64_C2(0xaa242499, 0x697392d3),
00187         RAPIDJSON_UINT64_C2(0xfd87b5f2, 0x8300ca0e), RAPIDJSON_UINT64_C2(0xbce50864, 0x92111aeb),
00188         RAPIDJSON_UINT64_C2(0x8cbccc09, 0x6f5088cc), RAPIDJSON_UINT64_C2(0xd1b71758, 0xe219652c),
00189         RAPIDJSON_UINT64_C2(0x9c400000, 0x00000000), RAPIDJSON_UINT64_C2(0xe8d4a510, 0x00000000),
00190         RAPIDJSON_UINT64_C2(0xad78ebc5, 0xac620000), RAPIDJSON_UINT64_C2(0x813f3978, 0xf8940984),
00191         RAPIDJSON_UINT64_C2(0xc097ce7b, 0xc90715b3), RAPIDJSON_UINT64_C2(0x8f7e32ce, 0x7bea5c70),
00192         RAPIDJSON_UINT64_C2(0xd5d238a4, 0xabe98068), RAPIDJSON_UINT64_C2(0x9f4f2726, 0x179a2245),
00193         RAPIDJSON_UINT64_C2(0xed63a231, 0xd4c4fb27), RAPIDJSON_UINT64_C2(0xb0de6538, 0x8cc8ada8),
00194         RAPIDJSON_UINT64_C2(0x83c7088e, 0x1aab65db), RAPIDJSON_UINT64_C2(0xc45d1df9, 0x42711d9a),
00195         RAPIDJSON_UINT64_C2(0x924d692c, 0xa61be758), RAPIDJSON_UINT64_C2(0xda01ee64, 0x1a708dea),
00196         RAPIDJSON_UINT64_C2(0xa26da399, 0x9aef774a), RAPIDJSON_UINT64_C2(0xf209787b, 0xb47d6b85),
00197         RAPIDJSON_UINT64_C2(0xb454e4a1, 0x79dd1877), RAPIDJSON_UINT64_C2(0x865b8692, 0x5b9bc5c2),
00198         RAPIDJSON_UINT64_C2(0xc83553c5, 0xc8965d3d), RAPIDJSON_UINT64_C2(0x952ab45c, 0xfa97a0b3),
00199         RAPIDJSON_UINT64_C2(0xde469fbd, 0x99a05fe3), RAPIDJSON_UINT64_C2(0xa59bc234, 0xdb398c25),
00200         RAPIDJSON_UINT64_C2(0xf6c69a72, 0xa3989f5c), RAPIDJSON_UINT64_C2(0xb7dcbf53, 0x54e9bece),
00201         RAPIDJSON_UINT64_C2(0x88fcf317, 0xf22241e2), RAPIDJSON_UINT64_C2(0xcc20ce9b, 0xd35c78a5),
00202         RAPIDJSON_UINT64_C2(0x98165af3, 0x7b2153df), RAPIDJSON_UINT64_C2(0xe2a0b5dc, 0x971f303a),
00203         RAPIDJSON_UINT64_C2(0xa8d9d153, 0x5ce3b396), RAPIDJSON_UINT64_C2(0xfb9b7cd9, 0xa4a7443c),
00204         RAPIDJSON_UINT64_C2(0xbb764c4c, 0xa7a44410), RAPIDJSON_UINT64_C2(0x8bab8eef, 0xb6409c1a),
00205         RAPIDJSON_UINT64_C2(0xd01fef10, 0xa657842c), RAPIDJSON_UINT64_C2(0x9b10a4e5, 0xe9913129),
00206         RAPIDJSON_UINT64_C2(0xe7109bfb, 0xa19c0c9d), RAPIDJSON_UINT64_C2(0xac2820d9, 0x623bf429),
00207         RAPIDJSON_UINT64_C2(0x80444b5e, 0x7aa7cf85), RAPIDJSON_UINT64_C2(0xbf21e440, 0x03acdd2d),
00208         RAPIDJSON_UINT64_C2(0x8e679c2f, 0x5e44ff8f), RAPIDJSON_UINT64_C2(0xd433179d, 0x9c8cb841),
00209         RAPIDJSON_UINT64_C2(0x9e19db92, 0xb4e31ba9), RAPIDJSON_UINT64_C2(0xeb96bf6e, 0xbadf77d9),
00210         RAPIDJSON_UINT64_C2(0xaf87023b, 0x9bf0ee6b)
00211     };
00212     static const int16_t kCachedPowers_E[] = {
00213         -1220, -1193, -1166, -1140, -1113, -1087, -1060, -1034, -1007,  -980,
00214         -954,  -927,  -901,  -874,  -847,  -821,  -794,  -768,  -741,  -715,
00215         -688,  -661,  -635,  -608,  -582,  -555,  -529,  -502,  -475,  -449,
00216         -422,  -396,  -369,  -343,  -316,  -289,  -263,  -236,  -210,  -183,
00217         -157,  -130,  -103,   -77,   -50,   -24,     3,    30,    56,    83,
00218         109,   136,   162,   189,   216,   242,   269,   295,   322,   348,
00219         375,   402,   428,   455,   481,   508,   534,   561,   588,   614,
00220         641,   667,   694,   720,   747,   774,   800,   827,   853,   880,
00221         907,   933,   960,   986,  1013,  1039,  1066
00222     };
00223     return DiyFp(kCachedPowers_F[index], kCachedPowers_E[index]);
00224 }
00225     
00226 inline DiyFp GetCachedPower(int e, int* K) {
00227 
00228     //int k = static_cast<int>(ceil((-61 - e) * 0.30102999566398114)) + 374;
00229     double dk = (-61 - e) * 0.30102999566398114 + 347;  // dk must be positive, so can do ceiling in positive
00230     int k = static_cast<int>(dk);
00231     if (dk - k > 0.0)
00232         k++;
00233 
00234     unsigned index = static_cast<unsigned>((k >> 3) + 1);
00235     *K = -(-348 + static_cast<int>(index << 3));    // decimal exponent no need lookup table
00236 
00237     return GetCachedPowerByIndex(index);
00238 }
00239 
00240 inline DiyFp GetCachedPower10(int exp, int *outExp) {
00241      unsigned index = (static_cast<unsigned>(exp) + 348u) / 8u;
00242      *outExp = -348 + static_cast<int>(index) * 8;
00243      return GetCachedPowerByIndex(index);
00244  }
00245 
00246 #ifdef __GNUC__
00247 RAPIDJSON_DIAG_POP
00248 #endif
00249 
00250 #ifdef __clang__
00251 RAPIDJSON_DIAG_POP
00252 RAPIDJSON_DIAG_OFF(padded)
00253 #endif
00254 
00255 } // namespace internal
00256 RAPIDJSON_NAMESPACE_END
00257 
00258 #endif // RAPIDJSON_DIYFP_H_