|
| 1 | +// Copyright 2016 Dtoa Developers |
| 2 | +// |
| 3 | +// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or |
| 4 | +// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license |
| 5 | +// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your |
| 6 | +// option. This file may not be copied, modified, or distributed |
| 7 | +// except according to those terms. |
| 8 | + |
| 9 | +use std::{mem, ops}; |
| 10 | + |
| 11 | +const DIY_SIGNIFICAND_SIZE: isize = 64; |
| 12 | +const DP_SIGNIFICAND_SIZE: isize = 52; |
| 13 | +const DP_EXPONENT_BIAS: isize = 0x3FF + DP_SIGNIFICAND_SIZE; |
| 14 | +const DP_MIN_EXPONENT: isize = -DP_EXPONENT_BIAS; |
| 15 | +const DP_EXPONENT_MASK: u64 = 0x7FF0000000000000; |
| 16 | +const DP_SIGNIFICAND_MASK: u64 = 0x000FFFFFFFFFFFFF; |
| 17 | +const DP_HIDDEN_BIT: u64 = 0x0010000000000000; |
| 18 | + |
| 19 | +#[derive(Copy, Clone)] |
| 20 | +pub struct DiyFp { |
| 21 | + pub f: u64, |
| 22 | + pub e: isize, |
| 23 | +} |
| 24 | + |
| 25 | +impl DiyFp { |
| 26 | + pub fn new(f: u64, e: isize) -> Self { |
| 27 | + DiyFp { f: f, e: e } |
| 28 | + } |
| 29 | + |
| 30 | + /* |
| 31 | + explicit DiyFp(double d) { |
| 32 | + union { |
| 33 | + double d; |
| 34 | + uint64_t u64; |
| 35 | + } u = { d }; |
| 36 | +
|
| 37 | + int biased_e = static_cast<int>((u.u64 & kDpExponentMask) >> kDpSignificandSize); |
| 38 | + uint64_t significand = (u.u64 & kDpSignificandMask); |
| 39 | + if (biased_e != 0) { |
| 40 | + f = significand + kDpHiddenBit; |
| 41 | + e = biased_e - kDpExponentBias; |
| 42 | + } |
| 43 | + else { |
| 44 | + f = significand; |
| 45 | + e = kDpMinExponent + 1; |
| 46 | + } |
| 47 | + } |
| 48 | + */ |
| 49 | + pub unsafe fn from_f64(d: f64) -> Self { |
| 50 | + let u: u64 = mem::transmute(d); |
| 51 | + |
| 52 | + let biased_e = ((u & DP_EXPONENT_MASK) >> DP_SIGNIFICAND_SIZE) as isize; |
| 53 | + let significand = u & DP_SIGNIFICAND_MASK; |
| 54 | + if biased_e != 0 { |
| 55 | + DiyFp { |
| 56 | + f: significand + DP_HIDDEN_BIT, |
| 57 | + e: biased_e - DP_EXPONENT_BIAS, |
| 58 | + } |
| 59 | + } else { |
| 60 | + DiyFp { |
| 61 | + f: significand, |
| 62 | + e: DP_MIN_EXPONENT + 1, |
| 63 | + } |
| 64 | + } |
| 65 | + } |
| 66 | + |
| 67 | + /* |
| 68 | + DiyFp Normalize() const { |
| 69 | + DiyFp res = *this; |
| 70 | + while (!(res.f & (static_cast<uint64_t>(1) << 63))) { |
| 71 | + res.f <<= 1; |
| 72 | + res.e--; |
| 73 | + } |
| 74 | + return res; |
| 75 | + } |
| 76 | + */ |
| 77 | + pub fn normalize(self) -> DiyFp { |
| 78 | + let mut res = self; |
| 79 | + while (res.f & (1u64 << 63)) == 0 { |
| 80 | + res.f <<= 1; |
| 81 | + res.e -= 1; |
| 82 | + } |
| 83 | + res |
| 84 | + } |
| 85 | + |
| 86 | + /* |
| 87 | + DiyFp NormalizeBoundary() const { |
| 88 | + DiyFp res = *this; |
| 89 | + while (!(res.f & (kDpHiddenBit << 1))) { |
| 90 | + res.f <<= 1; |
| 91 | + res.e--; |
| 92 | + } |
| 93 | + res.f <<= (kDiySignificandSize - kDpSignificandSize - 2); |
| 94 | + res.e = res.e - (kDiySignificandSize - kDpSignificandSize - 2); |
| 95 | + return res; |
| 96 | + } |
| 97 | + */ |
| 98 | + pub fn normalize_boundary(self) -> DiyFp { |
| 99 | + let mut res = self; |
| 100 | + while (res.f & DP_HIDDEN_BIT << 1) == 0 { |
| 101 | + res.f <<= 1; |
| 102 | + res.e -= 1; |
| 103 | + } |
| 104 | + res.f <<= DIY_SIGNIFICAND_SIZE - DP_SIGNIFICAND_SIZE - 2; |
| 105 | + res.e -= DIY_SIGNIFICAND_SIZE - DP_SIGNIFICAND_SIZE - 2; |
| 106 | + res |
| 107 | + } |
| 108 | + |
| 109 | + /* |
| 110 | + void NormalizedBoundaries(DiyFp* minus, DiyFp* plus) const { |
| 111 | + DiyFp pl = DiyFp((f << 1) + 1, e - 1).NormalizeBoundary(); |
| 112 | + DiyFp mi = (f == kDpHiddenBit) ? DiyFp((f << 2) - 1, e - 2) : DiyFp((f << 1) - 1, e - 1); |
| 113 | + mi.f <<= mi.e - pl.e; |
| 114 | + mi.e = pl.e; |
| 115 | + *plus = pl; |
| 116 | + *minus = mi; |
| 117 | + } |
| 118 | + */ |
| 119 | + pub fn normalized_boundaries(self) -> (DiyFp, DiyFp) { |
| 120 | + let pl = DiyFp::new((self.f << 1) + 1, self.e - 1).normalize_boundary(); |
| 121 | + let mut mi = if self.f == DP_HIDDEN_BIT { |
| 122 | + DiyFp::new((self.f << 2) - 1, self.e - 2) |
| 123 | + } else { |
| 124 | + DiyFp::new((self.f << 1) - 1, self.e - 1) |
| 125 | + }; |
| 126 | + mi.f <<= mi.e - pl.e; |
| 127 | + mi.e = pl.e; |
| 128 | + (mi, pl) |
| 129 | + } |
| 130 | +} |
| 131 | + |
| 132 | +impl ops::Sub for DiyFp { |
| 133 | + type Output = DiyFp; |
| 134 | + fn sub(self, rhs: DiyFp) -> DiyFp { |
| 135 | + DiyFp { |
| 136 | + f: self.f - rhs.f, |
| 137 | + e: self.e, |
| 138 | + } |
| 139 | + } |
| 140 | +} |
| 141 | + |
| 142 | +impl ops::Mul for DiyFp { |
| 143 | + type Output = DiyFp; |
| 144 | + fn mul(self, rhs: DiyFp) -> DiyFp { |
| 145 | + let m32 = 0xFFFFFFFFu64; |
| 146 | + let a = self.f >> 32; |
| 147 | + let b = self.f & m32; |
| 148 | + let c = rhs.f >> 32; |
| 149 | + let d = rhs.f & m32; |
| 150 | + let ac = a * c; |
| 151 | + let bc = b * c; |
| 152 | + let ad = a * d; |
| 153 | + let bd = b * d; |
| 154 | + let mut tmp = (bd >> 32) + (ad & m32) + (bc & m32); |
| 155 | + tmp += 1u64 << 31; // mult_round |
| 156 | + DiyFp { |
| 157 | + f: ac + (ad >> 32) + (bc >> 32) + (tmp >> 32), |
| 158 | + e: self.e + rhs.e + 64, |
| 159 | + } |
| 160 | + } |
| 161 | +} |
| 162 | + |
| 163 | +fn get_cached_power_by_index(index: usize) -> DiyFp { |
| 164 | + // 10^-348, 10^-340, ..., 10^340 |
| 165 | + static CACHED_POWERS_F: [u64; 87] = [ |
| 166 | + 0xfa8fd5a0081c0288, 0xbaaee17fa23ebf76, |
| 167 | + 0x8b16fb203055ac76, 0xcf42894a5dce35ea, |
| 168 | + 0x9a6bb0aa55653b2d, 0xe61acf033d1a45df, |
| 169 | + 0xab70fe17c79ac6ca, 0xff77b1fcbebcdc4f, |
| 170 | + 0xbe5691ef416bd60c, 0x8dd01fad907ffc3c, |
| 171 | + 0xd3515c2831559a83, 0x9d71ac8fada6c9b5, |
| 172 | + 0xea9c227723ee8bcb, 0xaecc49914078536d, |
| 173 | + 0x823c12795db6ce57, 0xc21094364dfb5637, |
| 174 | + 0x9096ea6f3848984f, 0xd77485cb25823ac7, |
| 175 | + 0xa086cfcd97bf97f4, 0xef340a98172aace5, |
| 176 | + 0xb23867fb2a35b28e, 0x84c8d4dfd2c63f3b, |
| 177 | + 0xc5dd44271ad3cdba, 0x936b9fcebb25c996, |
| 178 | + 0xdbac6c247d62a584, 0xa3ab66580d5fdaf6, |
| 179 | + 0xf3e2f893dec3f126, 0xb5b5ada8aaff80b8, |
| 180 | + 0x87625f056c7c4a8b, 0xc9bcff6034c13053, |
| 181 | + 0x964e858c91ba2655, 0xdff9772470297ebd, |
| 182 | + 0xa6dfbd9fb8e5b88f, 0xf8a95fcf88747d94, |
| 183 | + 0xb94470938fa89bcf, 0x8a08f0f8bf0f156b, |
| 184 | + 0xcdb02555653131b6, 0x993fe2c6d07b7fac, |
| 185 | + 0xe45c10c42a2b3b06, 0xaa242499697392d3, |
| 186 | + 0xfd87b5f28300ca0e, 0xbce5086492111aeb, |
| 187 | + 0x8cbccc096f5088cc, 0xd1b71758e219652c, |
| 188 | + 0x9c40000000000000, 0xe8d4a51000000000, |
| 189 | + 0xad78ebc5ac620000, 0x813f3978f8940984, |
| 190 | + 0xc097ce7bc90715b3, 0x8f7e32ce7bea5c70, |
| 191 | + 0xd5d238a4abe98068, 0x9f4f2726179a2245, |
| 192 | + 0xed63a231d4c4fb27, 0xb0de65388cc8ada8, |
| 193 | + 0x83c7088e1aab65db, 0xc45d1df942711d9a, |
| 194 | + 0x924d692ca61be758, 0xda01ee641a708dea, |
| 195 | + 0xa26da3999aef774a, 0xf209787bb47d6b85, |
| 196 | + 0xb454e4a179dd1877, 0x865b86925b9bc5c2, |
| 197 | + 0xc83553c5c8965d3d, 0x952ab45cfa97a0b3, |
| 198 | + 0xde469fbd99a05fe3, 0xa59bc234db398c25, |
| 199 | + 0xf6c69a72a3989f5c, 0xb7dcbf5354e9bece, |
| 200 | + 0x88fcf317f22241e2, 0xcc20ce9bd35c78a5, |
| 201 | + 0x98165af37b2153df, 0xe2a0b5dc971f303a, |
| 202 | + 0xa8d9d1535ce3b396, 0xfb9b7cd9a4a7443c, |
| 203 | + 0xbb764c4ca7a44410, 0x8bab8eefb6409c1a, |
| 204 | + 0xd01fef10a657842c, 0x9b10a4e5e9913129, |
| 205 | + 0xe7109bfba19c0c9d, 0xac2820d9623bf429, |
| 206 | + 0x80444b5e7aa7cf85, 0xbf21e44003acdd2d, |
| 207 | + 0x8e679c2f5e44ff8f, 0xd433179d9c8cb841, |
| 208 | + 0x9e19db92b4e31ba9, 0xeb96bf6ebadf77d9, |
| 209 | + 0xaf87023b9bf0ee6b, |
| 210 | + ]; |
| 211 | + static CACHED_POWERS_E: [i16; 87] = [ |
| 212 | + -1220, -1193, -1166, -1140, -1113, -1087, -1060, -1034, -1007, -980, |
| 213 | + -954, -927, -901, -874, -847, -821, -794, -768, -741, -715, |
| 214 | + -688, -661, -635, -608, -582, -555, -529, -502, -475, -449, |
| 215 | + -422, -396, -369, -343, -316, -289, -263, -236, -210, -183, |
| 216 | + -157, -130, -103, -77, -50, -24, 3, 30, 56, 83, |
| 217 | + 109, 136, 162, 189, 216, 242, 269, 295, 322, 348, |
| 218 | + 375, 402, 428, 455, 481, 508, 534, 561, 588, 614, |
| 219 | + 641, 667, 694, 720, 747, 774, 800, 827, 853, 880, |
| 220 | + 907, 933, 960, 986, 1013, 1039, 1066, |
| 221 | + ]; |
| 222 | + DiyFp::new(CACHED_POWERS_F[index], CACHED_POWERS_E[index] as isize) |
| 223 | +} |
| 224 | + |
| 225 | +/* |
| 226 | +inline DiyFp GetCachedPower(int e, int* K) { |
| 227 | + //int k = static_cast<int>(ceil((-61 - e) * 0.30102999566398114)) + 374; |
| 228 | + double dk = (-61 - e) * 0.30102999566398114 + 347; // dk must be positive, so can do ceiling in positive |
| 229 | + int k = static_cast<int>(dk); |
| 230 | + if (dk - k > 0.0) |
| 231 | + k++; |
| 232 | +
|
| 233 | + unsigned index = static_cast<unsigned>((k >> 3) + 1); |
| 234 | + *K = -(-348 + static_cast<int>(index << 3)); // decimal exponent no need lookup table |
| 235 | +
|
| 236 | + return GetCachedPowerByIndex(index); |
| 237 | +} |
| 238 | +*/ |
| 239 | +#[inline] |
| 240 | +pub fn get_cached_power(e: isize) -> (DiyFp, isize) { |
| 241 | + let dk = (-61 - e) as f64 * 0.30102999566398114f64 + 347f64; // dk must be positive, so can do ceiling in positive |
| 242 | + let mut k = dk as isize; |
| 243 | + if dk - k as f64 > 0.0 { |
| 244 | + k += 1; |
| 245 | + } |
| 246 | + |
| 247 | + let index = ((k >> 3) + 1) as usize; |
| 248 | + let k = -(-348 + (index << 3) as isize); // decimal exponent no need lookup table |
| 249 | + |
| 250 | + (get_cached_power_by_index(index), k) |
| 251 | +} |
0 commit comments