// Licensed under the Apache License, Version 2.0 or the MIT license // , at your // option. This file may not be copied, modified, or distributed // except according to those terms. // // --- // // The C++ implementation preserved here in comments is licensed as follows: // // Tencent is pleased to support the open source community by making RapidJSON // available. // // Copyright (C) 2015 THL A29 Limited, a Tencent company, and Milo Yip. All // rights reserved. // // Licensed under the MIT License (the "License"); you may not use this file // except in compliance with the License. You may obtain a copy of the License // at // // http://opensource.org/licenses/MIT // // Unless required by applicable law or agreed to in writing, software // distributed under the License is distributed on an "AS IS" BASIS, WITHOUT // WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the // License for the specific language governing permissions and limitations under // the License. use std::ops; #[derive(Copy, Clone, Debug)] pub struct DiyFp { pub f: F, pub e: E, } impl DiyFp { pub fn new(f: F, e: E) -> Self { DiyFp { f: f, e: e } } } impl ops::Mul for DiyFp { type Output = Self; fn mul(self, rhs: Self) -> Self { let mut tmp = self.f as u64 * rhs.f as u64; tmp += 1u64 << 31; // mult_round DiyFp { f: (tmp >> 32) as u32, e: self.e + rhs.e + 32, } } } impl ops::Mul for DiyFp { type Output = Self; fn mul(self, rhs: Self) -> Self { let m32 = 0xFFFFFFFFu64; let a = self.f >> 32; let b = self.f & m32; let c = rhs.f >> 32; let d = rhs.f & m32; let ac = a * c; let bc = b * c; let ad = a * d; let bd = b * d; let mut tmp = (bd >> 32) + (ad & m32) + (bc & m32); tmp += 1u64 << 31; // mult_round DiyFp { f: ac + (ad >> 32) + (bc >> 32) + (tmp >> 32), e: self.e + rhs.e + 64, } } } #[doc(hidden)] #[macro_export] macro_rules! diyfp {( floating_type: $fty:ty, significand_type: $sigty:ty, exponent_type: $expty:ty, diy_significand_size: $diy_significand_size:expr, significand_size: $significand_size:expr, exponent_bias: $exponent_bias:expr, mask_type: $mask_type:ty, exponent_mask: $exponent_mask:expr, significand_mask: $significand_mask:expr, hidden_bit: $hidden_bit:expr, cached_powers_f: $cached_powers_f:expr, cached_powers_e: $cached_powers_e:expr, min_power: $min_power:expr, ) => { type DiyFp = diyfp::DiyFp<$sigty, $expty>; impl DiyFp { // Preconditions: // `d` must have a positive sign and must not be infinity or NaN. /* explicit DiyFp(double d) { union { double d; uint64_t u64; } u = { d }; int biased_e = static_cast((u.u64 & kDpExponentMask) >> kDpSignificandSize); uint64_t significand = (u.u64 & kDpSignificandMask); if (biased_e != 0) { f = significand + kDpHiddenBit; e = biased_e - kDpExponentBias; } else { f = significand; e = kDpMinExponent + 1; } } */ unsafe fn from(d: $fty) -> Self { let u: $mask_type = mem::transmute(d); let biased_e = ((u & $exponent_mask) >> $significand_size) as $expty; let significand = u & $significand_mask; if biased_e != 0 { DiyFp { f: significand + $hidden_bit, e: biased_e - $exponent_bias - $significand_size, } } else { DiyFp { f: significand, e: 1 - $exponent_bias - $significand_size, } } } // Normalizes so that the highest bit of the diy significand is 1. /* DiyFp Normalize() const { DiyFp res = *this; while (!(res.f & (static_cast(1) << 63))) { res.f <<= 1; res.e--; } return res; } */ fn normalize(self) -> DiyFp { let mut res = self; while (res.f & (1 << ($diy_significand_size - 1))) == 0 { res.f <<= 1; res.e -= 1; } res } // Normalizes so that the highest bit of the diy significand is 1. // // Precondition: // `self.f` must be no more than 2 bits longer than the f64 significand. /* DiyFp NormalizeBoundary() const { DiyFp res = *this; while (!(res.f & (kDpHiddenBit << 1))) { res.f <<= 1; res.e--; } res.f <<= (kDiySignificandSize - kDpSignificandSize - 2); res.e = res.e - (kDiySignificandSize - kDpSignificandSize - 2); return res; } */ fn normalize_boundary(self) -> DiyFp { let mut res = self; while (res.f & $hidden_bit << 1) == 0 { res.f <<= 1; res.e -= 1; } res.f <<= $diy_significand_size - $significand_size - 2; res.e -= $diy_significand_size - $significand_size - 2; res } // Normalizes `self - e` and `self + e` where `e` is half of the least // significant digit of `self`. The plus is normalized so that the highest // bit of the diy significand is 1. The minus is normalized so that it has // the same exponent as the plus. // // Preconditions: // `self` must have been returned directly from `DiyFp::from_f64`. // `self.f` must not be zero. /* void NormalizedBoundaries(DiyFp* minus, DiyFp* plus) const { DiyFp pl = DiyFp((f << 1) + 1, e - 1).NormalizeBoundary(); DiyFp mi = (f == kDpHiddenBit) ? DiyFp((f << 2) - 1, e - 2) : DiyFp((f << 1) - 1, e - 1); mi.f <<= mi.e - pl.e; mi.e = pl.e; *plus = pl; *minus = mi; } */ fn normalized_boundaries(self) -> (DiyFp, DiyFp) { let pl = DiyFp::new((self.f << 1) + 1, self.e - 1).normalize_boundary(); let mut mi = if self.f == $hidden_bit { DiyFp::new((self.f << 2) - 1, self.e - 2) } else { DiyFp::new((self.f << 1) - 1, self.e - 1) }; mi.f <<= mi.e - pl.e; mi.e = pl.e; (mi, pl) } } impl ops::Sub for DiyFp { type Output = Self; fn sub(self, rhs: Self) -> Self { DiyFp { f: self.f - rhs.f, e: self.e, } } } /* inline DiyFp GetCachedPower(int e, int* K) { //int k = static_cast(ceil((-61 - e) * 0.30102999566398114)) + 374; double dk = (-61 - e) * 0.30102999566398114 + 347; // dk must be positive, so can do ceiling in positive int k = static_cast(dk); if (dk - k > 0.0) k++; unsigned index = static_cast((k >> 3) + 1); *K = -(-348 + static_cast(index << 3)); // decimal exponent no need lookup table return GetCachedPowerByIndex(index); } */ #[inline] fn get_cached_power(e: $expty) -> (DiyFp, isize) { let dk = (3 - $diy_significand_size - e) as f64 * 0.30102999566398114f64 - ($min_power + 1) as f64; let mut k = dk as isize; if dk - k as f64 > 0.0 { k += 1; } let index = ((k >> 3) + 1) as usize; let k = -($min_power + (index << 3) as isize); (DiyFp::new($cached_powers_f[index], $cached_powers_e[index] as $expty), k) } }}