1 // Copyright 2017 The Rust Project Developers. See the COPYRIGHT
2 // file at the top-level directory of this distribution and at
3 // http://rust-lang.org/COPYRIGHT.
5 // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
6 // http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
7 // <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
8 // option. This file may not be copied, modified, or distributed
9 // except according to those terms.
11 use {Category, ExpInt, Float, FloatConvert, Round, ParseError, Status, StatusAnd};
14 use std::cmp::Ordering;
19 #[derive(Copy, Clone, PartialEq, PartialOrd, Debug)]
20 pub struct DoubleFloat<F>(F, F);
21 pub type DoubleDouble = DoubleFloat<ieee::Double>;
23 // These are legacy semantics for the Fallback, inaccrurate implementation of
24 // IBM double-double, if the accurate DoubleDouble doesn't handle the
25 // operation. It's equivalent to having an IEEE number with consecutive 106
26 // bits of mantissa and 11 bits of exponent.
28 // It's not equivalent to IBM double-double. For example, a legit IBM
29 // double-double, 1 + epsilon:
31 // 1 + epsilon = 1 + (1 >> 1076)
33 // is not representable by a consecutive 106 bits of mantissa.
35 // Currently, these semantics are used in the following way:
37 // DoubleDouble -> (Double, Double) ->
38 // DoubleDouble's Fallback -> IEEE operations
40 // FIXME: Implement all operations in DoubleDouble, and delete these
42 // FIXME(eddyb) This shouldn't need to be `pub`, it's only used in bounds.
43 pub struct FallbackS<F>(F);
44 type Fallback<F> = ieee::IeeeFloat<FallbackS<F>>;
45 impl<F: Float> ieee::Semantics for FallbackS<F> {
46 // Forbid any conversion to/from bits.
47 const BITS: usize = 0;
48 const PRECISION: usize = F::PRECISION * 2;
49 const MAX_EXP: ExpInt = F::MAX_EXP as ExpInt;
50 const MIN_EXP: ExpInt = F::MIN_EXP as ExpInt + F::PRECISION as ExpInt;
53 // Convert number to F. To avoid spurious underflows, we re-
54 // normalize against the F exponent range first, and only *then*
55 // truncate the mantissa. The result of that second conversion
56 // may be inexact, but should never underflow.
57 // FIXME(eddyb) This shouldn't need to be `pub`, it's only used in bounds.
58 pub struct FallbackExtendedS<F>(F);
59 type FallbackExtended<F> = ieee::IeeeFloat<FallbackExtendedS<F>>;
60 impl<F: Float> ieee::Semantics for FallbackExtendedS<F> {
61 // Forbid any conversion to/from bits.
62 const BITS: usize = 0;
63 const PRECISION: usize = Fallback::<F>::PRECISION;
64 const MAX_EXP: ExpInt = F::MAX_EXP as ExpInt;
67 impl<F: Float> From<Fallback<F>> for DoubleFloat<F>
69 F: FloatConvert<FallbackExtended<F>>,
70 FallbackExtended<F>: FloatConvert<F>,
72 fn from(x: Fallback<F>) -> Self {
74 let mut loses_info = false;
76 let extended: FallbackExtended<F> = unpack!(status=, x.convert(&mut loses_info));
77 assert_eq!((status, loses_info), (Status::OK, false));
79 let a = unpack!(status=, extended.convert(&mut loses_info));
80 assert_eq!(status - Status::INEXACT, Status::OK);
82 // If conversion was exact or resulted in a special case, we're done;
83 // just set the second double to zero. Otherwise, re-convert back to
84 // the extended format and compute the difference. This now should
85 // convert exactly to double.
86 let b = if a.is_finite_non_zero() && loses_info {
87 let u: FallbackExtended<F> = unpack!(status=, a.convert(&mut loses_info));
88 assert_eq!((status, loses_info), (Status::OK, false));
89 let v = unpack!(status=, extended - u);
90 assert_eq!(status, Status::OK);
91 let v = unpack!(status=, v.convert(&mut loses_info));
92 assert_eq!((status, loses_info), (Status::OK, false));
102 impl<F: FloatConvert<Self>> From<DoubleFloat<F>> for Fallback<F> {
103 fn from(DoubleFloat(a, b): DoubleFloat<F>) -> Self {
105 let mut loses_info = false;
107 // Get the first F and convert to our format.
108 let a = unpack!(status=, a.convert(&mut loses_info));
109 assert_eq!((status, loses_info), (Status::OK, false));
111 // Unless we have a special case, add in second F.
112 if a.is_finite_non_zero() {
113 let b = unpack!(status=, b.convert(&mut loses_info));
114 assert_eq!((status, loses_info), (Status::OK, false));
123 float_common_impls!(DoubleFloat<F>);
125 impl<F: Float> Neg for DoubleFloat<F> {
127 fn neg(self) -> Self {
128 if self.1.is_finite_non_zero() {
129 DoubleFloat(-self.0, -self.1)
131 DoubleFloat(-self.0, self.1)
136 impl<F: FloatConvert<Fallback<F>>> fmt::Display for DoubleFloat<F> {
137 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
138 fmt::Display::fmt(&Fallback::from(*self), f)
142 impl<F: FloatConvert<Fallback<F>>> Float for DoubleFloat<F>
144 Self: From<Fallback<F>>,
146 const BITS: usize = F::BITS * 2;
147 const PRECISION: usize = Fallback::<F>::PRECISION;
148 const MAX_EXP: ExpInt = Fallback::<F>::MAX_EXP;
149 const MIN_EXP: ExpInt = Fallback::<F>::MIN_EXP;
151 const ZERO: Self = DoubleFloat(F::ZERO, F::ZERO);
153 const INFINITY: Self = DoubleFloat(F::INFINITY, F::ZERO);
155 // FIXME(eddyb) remove when qnan becomes const fn.
156 const NAN: Self = DoubleFloat(F::NAN, F::ZERO);
158 fn qnan(payload: Option<u128>) -> Self {
159 DoubleFloat(F::qnan(payload), F::ZERO)
162 fn snan(payload: Option<u128>) -> Self {
163 DoubleFloat(F::snan(payload), F::ZERO)
166 fn largest() -> Self {
168 let mut r = DoubleFloat(F::largest(), F::largest());
169 r.1 = r.1.scalbn(-(F::PRECISION as ExpInt + 1));
170 r.1 = unpack!(status=, r.1.next_down());
171 assert_eq!(status, Status::OK);
175 const SMALLEST: Self = DoubleFloat(F::SMALLEST, F::ZERO);
177 fn smallest_normalized() -> Self {
179 F::smallest_normalized().scalbn(F::PRECISION as ExpInt),
184 // Implement addition, subtraction, multiplication and division based on:
185 // "Software for Doubled-Precision Floating-Point Computations",
186 // by Seppo Linnainmaa, ACM TOMS vol 7 no 3, September 1981, pages 272-283.
188 fn add_r(mut self, rhs: Self, round: Round) -> StatusAnd<Self> {
189 match (self.category(), rhs.category()) {
190 (Category::Infinity, Category::Infinity) => {
191 if self.is_negative() != rhs.is_negative() {
192 Status::INVALID_OP.and(Self::NAN.copy_sign(self))
198 (_, Category::Zero) |
200 (Category::Infinity, Category::Normal) => Status::OK.and(self),
202 (Category::Zero, _) |
204 (_, Category::Infinity) => Status::OK.and(rhs),
206 (Category::Normal, Category::Normal) => {
207 let mut status = Status::OK;
208 let (a, aa, c, cc) = (self.0, self.1, rhs.0, rhs.1);
210 z = unpack!(status|=, z.add_r(c, round));
212 if !z.is_infinite() {
213 return status.and(DoubleFloat(z, F::ZERO));
216 let a_cmp_c = a.cmp_abs_normal(c);
218 z = unpack!(status|=, z.add_r(aa, round));
219 if a_cmp_c == Ordering::Greater {
220 // z = cc + aa + c + a;
221 z = unpack!(status|=, z.add_r(c, round));
222 z = unpack!(status|=, z.add_r(a, round));
224 // z = cc + aa + a + c;
225 z = unpack!(status|=, z.add_r(a, round));
226 z = unpack!(status|=, z.add_r(c, round));
229 return status.and(DoubleFloat(z, F::ZERO));
233 zz = unpack!(status|=, zz.add_r(cc, round));
234 if a_cmp_c == Ordering::Greater {
235 // self.1 = a - z + c + zz;
237 self.1 = unpack!(status|=, self.1.sub_r(z, round));
238 self.1 = unpack!(status|=, self.1.add_r(c, round));
239 self.1 = unpack!(status|=, self.1.add_r(zz, round));
241 // self.1 = c - z + a + zz;
243 self.1 = unpack!(status|=, self.1.sub_r(z, round));
244 self.1 = unpack!(status|=, self.1.add_r(a, round));
245 self.1 = unpack!(status|=, self.1.add_r(zz, round));
250 q = unpack!(status|=, q.sub_r(z, round));
252 // zz = q + c + (a - (q + z)) + aa + cc;
253 // Compute a - (q + z) as -((q + z) - a) to avoid temporary copies.
255 zz = unpack!(status|=, zz.add_r(c, round));
256 q = unpack!(status|=, q.add_r(z, round));
257 q = unpack!(status|=, q.sub_r(a, round));
259 zz = unpack!(status|=, zz.add_r(q, round));
260 zz = unpack!(status|=, zz.add_r(aa, round));
261 zz = unpack!(status|=, zz.add_r(cc, round));
262 if zz.is_zero() && !zz.is_negative() {
263 return Status::OK.and(DoubleFloat(z, F::ZERO));
266 self.0 = unpack!(status|=, self.0.add_r(zz, round));
267 if !self.0.is_finite() {
269 return status.and(self);
272 self.1 = unpack!(status|=, self.1.sub_r(self.0, round));
273 self.1 = unpack!(status|=, self.1.add_r(zz, round));
280 fn mul_r(mut self, rhs: Self, round: Round) -> StatusAnd<Self> {
281 // Interesting observation: For special categories, finding the lowest
282 // common ancestor of the following layered graph gives the correct
291 // e.g. NaN * NaN = NaN
293 // Normal * Zero = Zero
294 // Normal * Inf = Inf
295 match (self.category(), rhs.category()) {
296 (Category::NaN, _) => Status::OK.and(self),
298 (_, Category::NaN) => Status::OK.and(rhs),
300 (Category::Zero, Category::Infinity) |
301 (Category::Infinity, Category::Zero) => Status::OK.and(Self::NAN),
303 (Category::Zero, _) |
304 (Category::Infinity, _) => Status::OK.and(self),
306 (_, Category::Zero) |
307 (_, Category::Infinity) => Status::OK.and(rhs),
309 (Category::Normal, Category::Normal) => {
310 let mut status = Status::OK;
311 let (a, b, c, d) = (self.0, self.1, rhs.0, rhs.1);
314 t = unpack!(status|=, t.mul_r(c, round));
315 if !t.is_finite_non_zero() {
316 return status.and(DoubleFloat(t, F::ZERO));
319 // tau = fmsub(a, c, t), that is -fmadd(-a, c, t).
321 tau = unpack!(status|=, tau.mul_add_r(c, -t, round));
324 v = unpack!(status|=, v.mul_r(d, round));
327 w = unpack!(status|=, w.mul_r(c, round));
328 v = unpack!(status|=, v.add_r(w, round));
330 tau = unpack!(status|=, tau.add_r(v, round));
333 u = unpack!(status|=, u.add_r(tau, round));
339 // self.1 = (t - u) + tau
340 t = unpack!(status|=, t.sub_r(u, round));
341 t = unpack!(status|=, t.add_r(tau, round));
349 fn mul_add_r(self, multiplicand: Self, addend: Self, round: Round) -> StatusAnd<Self> {
351 .mul_add_r(Fallback::from(multiplicand), Fallback::from(addend), round)
355 fn div_r(self, rhs: Self, round: Round) -> StatusAnd<Self> {
356 Fallback::from(self).div_r(Fallback::from(rhs), round).map(
361 fn c_fmod(self, rhs: Self) -> StatusAnd<Self> {
362 Fallback::from(self).c_fmod(Fallback::from(rhs)).map(
367 fn round_to_integral(self, round: Round) -> StatusAnd<Self> {
368 Fallback::from(self).round_to_integral(round).map(
373 fn next_up(self) -> StatusAnd<Self> {
374 Fallback::from(self).next_up().map(Self::from)
377 fn from_bits(input: u128) -> Self {
378 let (a, b) = (input, input >> F::BITS);
380 F::from_bits(a & ((1 << F::BITS) - 1)),
381 F::from_bits(b & ((1 << F::BITS) - 1)),
385 fn from_u128_r(input: u128, round: Round) -> StatusAnd<Self> {
386 Fallback::from_u128_r(input, round).map(Self::from)
389 fn from_str_r(s: &str, round: Round) -> Result<StatusAnd<Self>, ParseError> {
390 Fallback::from_str_r(s, round).map(|r| r.map(Self::from))
393 fn to_bits(self) -> u128 {
394 self.0.to_bits() | (self.1.to_bits() << F::BITS)
397 fn to_u128_r(self, width: usize, round: Round, is_exact: &mut bool) -> StatusAnd<u128> {
398 Fallback::from(self).to_u128_r(width, round, is_exact)
401 fn cmp_abs_normal(self, rhs: Self) -> Ordering {
402 self.0.cmp_abs_normal(rhs.0).then_with(|| {
403 let result = self.1.cmp_abs_normal(rhs.1);
404 if result != Ordering::Equal {
405 let against = self.0.is_negative() ^ self.1.is_negative();
406 let rhs_against = rhs.0.is_negative() ^ rhs.1.is_negative();
407 (!against).cmp(&!rhs_against).then_with(|| if against {
418 fn bitwise_eq(self, rhs: Self) -> bool {
419 self.0.bitwise_eq(rhs.0) && self.1.bitwise_eq(rhs.1)
422 fn is_negative(self) -> bool {
426 fn is_denormal(self) -> bool {
427 self.category() == Category::Normal &&
428 (self.0.is_denormal() || self.0.is_denormal() ||
429 // (double)(Hi + Lo) == Hi defines a normal number.
430 !(self.0 + self.1).value.bitwise_eq(self.0))
433 fn is_signaling(self) -> bool {
434 self.0.is_signaling()
437 fn category(self) -> Category {
441 fn get_exact_inverse(self) -> Option<Self> {
442 Fallback::from(self).get_exact_inverse().map(Self::from)
445 fn ilogb(self) -> ExpInt {
449 fn scalbn_r(self, exp: ExpInt, round: Round) -> Self {
450 DoubleFloat(self.0.scalbn_r(exp, round), self.1.scalbn_r(exp, round))
453 fn frexp_r(self, exp: &mut ExpInt, round: Round) -> Self {
454 let a = self.0.frexp_r(exp, round);
456 if self.category() == Category::Normal {
457 b = b.scalbn_r(-*exp, round);