1 // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or 2 // http://www.apache.org/licenses/LICENSE-2.0> or the MIT license 3 // <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your 4 // option. This file may not be copied, modified, or distributed 5 // except according to those terms. 6 // 7 // --- 8 // 9 // The C++ implementation preserved here in comments is licensed as follows: 10 // 11 // Tencent is pleased to support the open source community by making RapidJSON 12 // available. 13 // 14 // Copyright (C) 2015 THL A29 Limited, a Tencent company, and Milo Yip. All 15 // rights reserved. 16 // 17 // Licensed under the MIT License (the "License"); you may not use this file 18 // except in compliance with the License. You may obtain a copy of the License 19 // at 20 // 21 // http://opensource.org/licenses/MIT 22 // 23 // Unless required by applicable law or agreed to in writing, software 24 // distributed under the License is distributed on an "AS IS" BASIS, WITHOUT 25 // WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the 26 // License for the specific language governing permissions and limitations under 27 // the License. 28 29 use std::ops; 30 31 #[derive(Copy, Clone, Debug)] 32 pub struct DiyFp<F, E> { 33 pub f: F, 34 pub e: E, 35 } 36 37 impl<F, E> DiyFp<F, E> { new(f: F, e: E) -> Self38 pub fn new(f: F, e: E) -> Self { 39 DiyFp { f: f, e: e } 40 } 41 } 42 43 impl ops::Mul for DiyFp<u32, i32> { 44 type Output = Self; mul(self, rhs: Self) -> Self45 fn mul(self, rhs: Self) -> Self { 46 let mut tmp = self.f as u64 * rhs.f as u64; 47 tmp += 1u64 << 31; // mult_round 48 DiyFp { 49 f: (tmp >> 32) as u32, 50 e: self.e + rhs.e + 32, 51 } 52 } 53 } 54 55 impl ops::Mul for DiyFp<u64, isize> { 56 type Output = Self; mul(self, rhs: Self) -> Self57 fn mul(self, rhs: Self) -> Self { 58 let m32 = 0xFFFFFFFFu64; 59 let a = self.f >> 32; 60 let b = self.f & m32; 61 let c = rhs.f >> 32; 62 let d = rhs.f & m32; 63 let ac = a * c; 64 let bc = b * c; 65 let ad = a * d; 66 let bd = b * d; 67 let mut tmp = (bd >> 32) + (ad & m32) + (bc & m32); 68 tmp += 1u64 << 31; // mult_round 69 DiyFp { 70 f: ac + (ad >> 32) + (bc >> 32) + (tmp >> 32), 71 e: self.e + rhs.e + 64, 72 } 73 } 74 } 75 76 #[doc(hidden)] 77 #[macro_export] 78 macro_rules! diyfp {( 79 floating_type: $fty:ty, 80 significand_type: $sigty:ty, 81 exponent_type: $expty:ty, 82 83 diy_significand_size: $diy_significand_size:expr, 84 significand_size: $significand_size:expr, 85 exponent_bias: $exponent_bias:expr, 86 mask_type: $mask_type:ty, 87 exponent_mask: $exponent_mask:expr, 88 significand_mask: $significand_mask:expr, 89 hidden_bit: $hidden_bit:expr, 90 cached_powers_f: $cached_powers_f:expr, 91 cached_powers_e: $cached_powers_e:expr, 92 min_power: $min_power:expr, 93 ) => { 94 95 type DiyFp = diyfp::DiyFp<$sigty, $expty>; 96 97 impl DiyFp { 98 // Preconditions: 99 // `d` must have a positive sign and must not be infinity or NaN. 100 /* 101 explicit DiyFp(double d) { 102 union { 103 double d; 104 uint64_t u64; 105 } u = { d }; 106 107 int biased_e = static_cast<int>((u.u64 & kDpExponentMask) >> kDpSignificandSize); 108 uint64_t significand = (u.u64 & kDpSignificandMask); 109 if (biased_e != 0) { 110 f = significand + kDpHiddenBit; 111 e = biased_e - kDpExponentBias; 112 } 113 else { 114 f = significand; 115 e = kDpMinExponent + 1; 116 } 117 } 118 */ 119 unsafe fn from(d: $fty) -> Self { 120 let u: $mask_type = mem::transmute(d); 121 122 let biased_e = ((u & $exponent_mask) >> $significand_size) as $expty; 123 let significand = u & $significand_mask; 124 if biased_e != 0 { 125 DiyFp { 126 f: significand + $hidden_bit, 127 e: biased_e - $exponent_bias - $significand_size, 128 } 129 } else { 130 DiyFp { 131 f: significand, 132 e: 1 - $exponent_bias - $significand_size, 133 } 134 } 135 } 136 137 // Normalizes so that the highest bit of the diy significand is 1. 138 /* 139 DiyFp Normalize() const { 140 DiyFp res = *this; 141 while (!(res.f & (static_cast<uint64_t>(1) << 63))) { 142 res.f <<= 1; 143 res.e--; 144 } 145 return res; 146 } 147 */ 148 fn normalize(self) -> DiyFp { 149 let mut res = self; 150 while (res.f & (1 << ($diy_significand_size - 1))) == 0 { 151 res.f <<= 1; 152 res.e -= 1; 153 } 154 res 155 } 156 157 // Normalizes so that the highest bit of the diy significand is 1. 158 // 159 // Precondition: 160 // `self.f` must be no more than 2 bits longer than the f64 significand. 161 /* 162 DiyFp NormalizeBoundary() const { 163 DiyFp res = *this; 164 while (!(res.f & (kDpHiddenBit << 1))) { 165 res.f <<= 1; 166 res.e--; 167 } 168 res.f <<= (kDiySignificandSize - kDpSignificandSize - 2); 169 res.e = res.e - (kDiySignificandSize - kDpSignificandSize - 2); 170 return res; 171 } 172 */ 173 fn normalize_boundary(self) -> DiyFp { 174 let mut res = self; 175 while (res.f & $hidden_bit << 1) == 0 { 176 res.f <<= 1; 177 res.e -= 1; 178 } 179 res.f <<= $diy_significand_size - $significand_size - 2; 180 res.e -= $diy_significand_size - $significand_size - 2; 181 res 182 } 183 184 // Normalizes `self - e` and `self + e` where `e` is half of the least 185 // significant digit of `self`. The plus is normalized so that the highest 186 // bit of the diy significand is 1. The minus is normalized so that it has 187 // the same exponent as the plus. 188 // 189 // Preconditions: 190 // `self` must have been returned directly from `DiyFp::from_f64`. 191 // `self.f` must not be zero. 192 /* 193 void NormalizedBoundaries(DiyFp* minus, DiyFp* plus) const { 194 DiyFp pl = DiyFp((f << 1) + 1, e - 1).NormalizeBoundary(); 195 DiyFp mi = (f == kDpHiddenBit) ? DiyFp((f << 2) - 1, e - 2) : DiyFp((f << 1) - 1, e - 1); 196 mi.f <<= mi.e - pl.e; 197 mi.e = pl.e; 198 *plus = pl; 199 *minus = mi; 200 } 201 */ 202 fn normalized_boundaries(self) -> (DiyFp, DiyFp) { 203 let pl = DiyFp::new((self.f << 1) + 1, self.e - 1).normalize_boundary(); 204 let mut mi = if self.f == $hidden_bit { 205 DiyFp::new((self.f << 2) - 1, self.e - 2) 206 } else { 207 DiyFp::new((self.f << 1) - 1, self.e - 1) 208 }; 209 mi.f <<= mi.e - pl.e; 210 mi.e = pl.e; 211 (mi, pl) 212 } 213 } 214 215 impl ops::Sub for DiyFp { 216 type Output = Self; 217 fn sub(self, rhs: Self) -> Self { 218 DiyFp { 219 f: self.f - rhs.f, 220 e: self.e, 221 } 222 } 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 fn get_cached_power(e: $expty) -> (DiyFp, isize) { 241 let dk = (3 - $diy_significand_size - e) as f64 * 0.30102999566398114f64 - ($min_power + 1) as f64; 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 = -($min_power + (index << 3) as isize); 249 250 (DiyFp::new($cached_powers_f[index], $cached_powers_e[index] as $expty), k) 251 } 252 253 }} 254