competitive/num/
double_double.rs

1use super::{Bounded, Decimal, IterScan, One, Zero};
2use std::{
3    cmp::Ordering,
4    fmt::{self, Display},
5    num::ParseFloatError,
6    ops::{Add, Div, Mul, Neg, Sub},
7    str::FromStr,
8};
9
10#[derive(Clone, Copy, Debug, Default, PartialEq)]
11pub struct DoubleDouble(f64, f64);
12
13impl DoubleDouble {
14    fn renormalize(a0: f64, a1: f64, a2: f64) -> Self {
15        let (s, t2) = quick_two_sum(a1, a2);
16        let (mut s, t1) = quick_two_sum(a0, s);
17        let mut b = (s, t1);
18        if t1 != 0. {
19            b.0 = s;
20            s = t1;
21        }
22        let (s, e) = quick_two_sum(s, t2);
23        if e != 0. {
24            if t1 != 0. {
25                b.1 = s;
26            } else {
27                b.0 = s;
28            }
29        }
30        Self(b.0, b.1)
31    }
32}
33
34fn quick_two_sum(a: f64, b: f64) -> (f64, f64) {
35    let s = a + b;
36    let e = b - (s - a);
37    (s, e)
38}
39
40fn two_sum(a: f64, b: f64) -> (f64, f64) {
41    let s = a + b;
42    let v = s - a;
43    let e = (a - (s - v)) + (b - v);
44    (s, e)
45}
46
47fn split(a: f64) -> (f64, f64) {
48    let t = 134_217_729. * a; // 134217729 = 2 ** 27 + 1
49    let ahi = t - (t - a);
50    let alo = a - ahi;
51    (ahi, alo)
52}
53
54fn two_prod(a: f64, b: f64) -> (f64, f64) {
55    let p = a * b;
56    let (ahi, alo) = split(a);
57    let (bhi, blo) = split(b);
58    let e = ((ahi * bhi - p) + ahi * blo + alo * bhi) + alo * blo;
59    (p, e)
60}
61
62fn three_two_sum(a: f64, b: f64, c: f64) -> (f64, f64) {
63    let (u, v) = two_sum(a, b);
64    let (r0, w) = two_sum(u, c);
65    let r1 = v + w;
66    (r0, r1)
67}
68
69impl Add<f64> for DoubleDouble {
70    type Output = Self;
71    fn add(self, rhs: f64) -> Self::Output {
72        let (t0, e) = two_sum(self.0, rhs);
73        let (t1, t2) = two_sum(self.1, e);
74        Self::renormalize(t0, t1, t2)
75    }
76}
77
78impl Add<DoubleDouble> for DoubleDouble {
79    type Output = Self;
80    fn add(self, rhs: Self) -> Self::Output {
81        let (t0, e) = two_sum(self.0, rhs.0);
82        let (t1, t2) = three_two_sum(self.1, rhs.1, e);
83        Self::renormalize(t0, t1, t2)
84    }
85}
86
87impl Sub for DoubleDouble {
88    type Output = Self;
89    fn sub(self, rhs: Self) -> Self::Output {
90        self + -rhs
91    }
92}
93
94impl Neg for DoubleDouble {
95    type Output = Self;
96    fn neg(self) -> Self::Output {
97        Self(-self.0, -self.1)
98    }
99}
100
101impl Mul<f64> for DoubleDouble {
102    type Output = Self;
103    fn mul(self, rhs: f64) -> Self::Output {
104        let (t0, e0) = two_prod(self.0, rhs);
105        let p1 = self.1 * rhs;
106        let (t1, t2) = two_sum(p1, e0);
107        Self::renormalize(t0, t1, t2)
108    }
109}
110
111impl Mul<DoubleDouble> for DoubleDouble {
112    type Output = Self;
113    fn mul(self, rhs: Self) -> Self::Output {
114        let (t0, q00) = two_prod(self.0, rhs.0);
115        let (p01, q01) = two_prod(self.0, rhs.1);
116        let (p10, q10) = two_prod(self.1, rhs.0);
117        let p11 = self.1 * rhs.1;
118        let (t1, e1) = three_two_sum(q00, p01, p10);
119        let t2 = e1 + q01 + q10 + p11;
120        Self::renormalize(t0, t1, t2)
121    }
122}
123
124impl Div<DoubleDouble> for DoubleDouble {
125    type Output = Self;
126    fn div(self, rhs: Self) -> Self::Output {
127        let q0 = self.0 / rhs.0;
128        let r = self - rhs * q0;
129        let q1 = r.0 / rhs.0;
130        let r = r - rhs * q1;
131        let q2 = r.0 / rhs.0;
132        Self::renormalize(q0, q1, q2)
133    }
134}
135
136impl From<DoubleDouble> for f64 {
137    fn from(x: DoubleDouble) -> f64 {
138        x.1 + x.0
139    }
140}
141
142impl From<DoubleDouble> for i64 {
143    fn from(mut x: DoubleDouble) -> i64 {
144        let is_neg = x.0.is_sign_negative();
145        if is_neg {
146            x = -x;
147        }
148        let mut i = 0i64;
149        for k in (1..64).rev() {
150            let t = (k as f64).exp2();
151            if x.0 >= t {
152                x = x + -t;
153                i += 1 << k;
154            }
155        }
156        i += x.0.round() as i64;
157        if is_neg {
158            i = -i;
159        }
160        i
161    }
162}
163
164impl From<f64> for DoubleDouble {
165    fn from(x: f64) -> Self {
166        Self(x, 0.)
167    }
168}
169
170impl Display for DoubleDouble {
171    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
172        write!(f, "{}", Decimal::from(self.0) + Decimal::from(self.1))
173    }
174}
175
176#[derive(Debug, Clone)]
177pub enum ParseDoubleDoubleError {
178    ParseFloatError(ParseFloatError),
179    ParseDecimalError(super::decimal::convert::ParseDecimalError),
180}
181
182impl From<ParseFloatError> for ParseDoubleDoubleError {
183    fn from(e: ParseFloatError) -> Self {
184        Self::ParseFloatError(e)
185    }
186}
187
188impl From<super::decimal::convert::ParseDecimalError> for ParseDoubleDoubleError {
189    fn from(e: super::decimal::convert::ParseDecimalError) -> Self {
190        Self::ParseDecimalError(e)
191    }
192}
193
194impl FromStr for DoubleDouble {
195    type Err = ParseDoubleDoubleError;
196    fn from_str(s: &str) -> Result<Self, Self::Err> {
197        let f0: f64 = s.parse()?;
198        let d1 = Decimal::from_str(s)? - Decimal::from(f0);
199        let f1: f64 = d1.to_string().parse()?;
200        Ok(Self::renormalize(f0, f1, 0.))
201    }
202}
203
204impl Eq for DoubleDouble {}
205impl PartialOrd for DoubleDouble {
206    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
207        Some(self.cmp(other))
208    }
209}
210impl Ord for DoubleDouble {
211    fn cmp(&self, other: &Self) -> Ordering {
212        fn total_cmp(x: f64, y: f64) -> Ordering {
213            let mut left = x.to_bits() as i64;
214            let mut right = y.to_bits() as i64;
215            left ^= (((left >> 63) as u64) >> 1) as i64;
216            right ^= (((right >> 63) as u64) >> 1) as i64;
217            left.cmp(&right)
218        }
219        total_cmp(self.0, other.0).then_with(|| total_cmp(self.1, other.1))
220    }
221}
222impl Bounded for DoubleDouble {
223    fn maximum() -> Self {
224        DoubleDouble::from(<f64 as Bounded>::maximum())
225    }
226    fn minimum() -> Self {
227        DoubleDouble::from(<f64 as Bounded>::minimum())
228    }
229}
230
231impl Zero for DoubleDouble {
232    fn zero() -> Self {
233        Self::from(0.)
234    }
235    fn is_zero(&self) -> bool
236    where
237        Self: PartialEq,
238    {
239        self.0 == 0.
240    }
241}
242
243impl One for DoubleDouble {
244    fn one() -> Self {
245        Self::from(1.)
246    }
247    fn is_one(&self) -> bool
248    where
249        Self: PartialEq,
250    {
251        self.0 == 1.
252    }
253}
254
255impl IterScan for DoubleDouble {
256    type Output = Self;
257    fn scan<'a, I: Iterator<Item = &'a str>>(iter: &mut I) -> Option<Self::Output> {
258        iter.next().and_then(|s| s.parse().ok())
259    }
260}
261
262impl DoubleDouble {
263    pub fn sqrt(self) -> Self {
264        if self.is_zero() {
265            return Self::from(0.);
266        }
267        let x = Self::from(1. / self.0.sqrt());
268        let x = x + x * (Self::from(1.) - self * x * x).div2(2.);
269        let x = x + x * (Self::from(1.) - self * x * x).div2(2.);
270        let x = x + x * (Self::from(1.) - self * x * x).div2(2.);
271        x * self
272    }
273    pub fn abs(self) -> Self {
274        if self.0.is_sign_negative() {
275            -self
276        } else {
277            self
278        }
279    }
280    fn div2(self, rhs: f64) -> Self {
281        Self(self.0 / rhs, self.1 / rhs)
282    }
283}
284
285#[cfg(test)]
286mod tests {
287    use super::*;
288
289    #[test]
290    fn test_display() {
291        let x = DoubleDouble::from(1.234);
292        assert_eq!(x.to_string(), "1.234");
293        let x = DoubleDouble::from(1.234e-10);
294        assert_eq!(x.to_string(), "0.0000000001234");
295        let x = DoubleDouble::from(1.234e10);
296        assert_eq!(x.to_string(), "12340000000");
297        let x = DoubleDouble::from(1.234e-10) + DoubleDouble::from(1.234e10);
298        assert_eq!(x.to_string(), "12340000000.0000000001234");
299    }
300
301    #[test]
302    fn test_from_str() {
303        let x = DoubleDouble::from_str("1.234").unwrap();
304        assert_eq!(x, DoubleDouble::from(1.234));
305        let x = DoubleDouble::from_str("0.0000000001234").unwrap();
306        assert_eq!(x, DoubleDouble::from(1.234e-10));
307        let x = DoubleDouble::from_str("12340000000").unwrap();
308        assert_eq!(x, DoubleDouble::from(1.234e10));
309        let x = DoubleDouble::from_str("12340000000.0000000001234").unwrap();
310        assert_eq!(
311            x,
312            DoubleDouble::from(1.234e10) + DoubleDouble::from(1.234e-10)
313        );
314    }
315}
competitive/num/double_double.rs

competitive/num/
double_double.rs
