competitive/num/
double_double.rs

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

competitive/num/
double_double.rs
