Struct QuadDouble

Source

pub struct QuadDouble(f64, f64, f64, f64);

Expand description

ref: https://na-inet.jp/na/qd_ja.pdf

Tuple Fields§

§0: f64§1: f64§2: f64§3: f64

Implementations§

Source §

impl QuadDouble

Source

fn renormalize(a0: f64, a1: f64, a2: f64, a3: f64, a4: f64) -> Self

Examples found in repository ?

crates/competitive/src/num/quad_double.rs (line 106)

101    fn add(self, rhs: f64) -> Self::Output {
102        let (t0, e) = two_sum(self.0, rhs);
103        let (t1, e) = two_sum(self.1, e);
104        let (t2, e) = two_sum(self.2, e);
105        let (t3, t4) = two_sum(self.3, e);
106        Self::renormalize(t0, t1, t2, t3, t4)
107    }
108}
109
110fn double_accumulate(u: f64, v: f64, x: f64) -> (f64, f64, f64) {
111    let (s, mut v) = two_sum(v, x);
112    let (mut s, mut u) = two_sum(u, s);
113    if u == 0. {
114        u = s;
115        s = 0.;
116    }
117    if v == 0. {
118        v = u;
119        u = s;
120        s = 0.
121    }
122    (s, u, v)
123}
124
125impl Add<QuadDouble> for QuadDouble {
126    type Output = Self;
127    fn add(self, rhs: Self) -> Self::Output {
128        let mut x = [0.; 8];
129        let (mut i, mut j, mut k) = (0, 0, 0);
130        while k < 8 {
131            if j >= 4 || i < 4 && self[i].abs() > rhs[j].abs() {
132                x[k] = self[i];
133                i += 1;
134            } else {
135                x[k] = rhs[j];
136                j += 1;
137            }
138            k += 1;
139        }
140
141        let (mut u, mut v) = (0., 0.);
142        let (mut k, mut i) = (0, 0);
143        let mut c = [0.; 4];
144        while k < 4 && i < 8 {
145            let tpl = double_accumulate(u, v, x[i]);
146            let s = tpl.0;
147            u = tpl.1;
148            v = tpl.2;
149            if s != 0. {
150                c[k] = s;
151                k += 1;
152            }
153            i += 1;
154        }
155        if k < 2 {
156            c[k + 1] = v;
157        }
158        if k < 3 {
159            c[k] = u;
160        }
161        Self::renormalize(c[0], c[1], c[2], c[3], 0.)
162    }
163}
164
165impl Sub for QuadDouble {
166    type Output = Self;
167    fn sub(self, rhs: Self) -> Self::Output {
168        self + -rhs
169    }
170}
171
172impl Neg for QuadDouble {
173    type Output = Self;
174    fn neg(self) -> Self::Output {
175        Self(-self.0, -self.1, -self.2, -self.3)
176    }
177}
178
179impl Mul<f64> for QuadDouble {
180    type Output = Self;
181    fn mul(self, rhs: f64) -> Self::Output {
182        let (t0, e0) = two_prod(self.0, rhs);
183        let (p1, e1) = two_prod(self.1, rhs);
184        let (p2, e2) = two_prod(self.2, rhs);
185        let p3 = self.3 * rhs;
186
187        let (t1, e4) = two_sum(p1, e0);
188        let (t2, e5, e6) = three_three_sum(p2, e1, e4);
189        let (t3, e7) = three_two_sum(p3, e2, e5);
190        let t4 = e7 + e6;
191        Self::renormalize(t0, t1, t2, t3, t4)
192    }
193}
194
195impl Mul<QuadDouble> for QuadDouble {
196    type Output = Self;
197    fn mul(self, rhs: Self) -> Self::Output {
198        let (t0, q00) = two_prod(self.0, rhs.0);
199
200        let (p01, q01) = two_prod(self.0, rhs.1);
201        let (p10, q10) = two_prod(self.1, rhs.0);
202
203        let (p02, q02) = two_prod(self.0, rhs.2);
204        let (p11, q11) = two_prod(self.1, rhs.1);
205        let (p20, q20) = two_prod(self.2, rhs.0);
206
207        let (p03, q03) = two_prod(self.0, rhs.3);
208        let (p12, q12) = two_prod(self.1, rhs.2);
209        let (p21, q21) = two_prod(self.2, rhs.1);
210        let (p30, q30) = two_prod(self.3, rhs.0);
211
212        let p13 = self.1 * rhs.3;
213        let p22 = self.2 * rhs.2;
214        let p31 = self.3 * rhs.1;
215
216        let (t1, e1, e2) = three_three_sum(q00, p01, p10);
217        let (t2, e3, e4) = multiple_three_sum(&[e1, q01, q10, p02, p11, p20]);
218        let (t3, e5) = multiple_two_sum(&[e2, e3, q02, q11, q20, p03, p12, p21, p30]);
219        let t4 = e4 + e5 + q03 + q12 + q21 + q30 + p13 + p22 + p31;
220        Self::renormalize(t0, t1, t2, t3, t4)
221    }
222}
223
224impl Div<QuadDouble> for QuadDouble {
225    type Output = Self;
226    fn div(self, rhs: Self) -> Self::Output {
227        let q0 = self.0 / rhs.0;
228        let r = self - rhs * q0;
229        let q1 = r.0 / rhs.0;
230        let r = r - rhs * q1;
231        let q2 = r.0 / rhs.0;
232        let r = r - rhs * q2;
233        let q3 = r.0 / rhs.0;
234        let r = r - rhs * q3;
235        let q4 = r.0 / rhs.0;
236        Self::renormalize(q0, q1, q2, q3, q4)
237    }
238}
239
240impl Index<usize> for QuadDouble {
241    type Output = f64;
242    fn index(&self, index: usize) -> &Self::Output {
243        match index {
244            0 => &self.0,
245            1 => &self.1,
246            2 => &self.2,
247            3 => &self.3,
248            _ => panic!(),
249        }
250    }
251}
252
253impl From<QuadDouble> for f64 {
254    fn from(x: QuadDouble) -> f64 {
255        x.3 + x.2 + x.1 + x.0
256    }
257}
258
259impl From<QuadDouble> for i64 {
260    fn from(mut x: QuadDouble) -> i64 {
261        let is_neg = x.0.is_sign_negative();
262        if is_neg {
263            x = -x;
264        }
265        let mut i = 0i64;
266        for k in (1..64).rev() {
267            let t = (k as f64).exp2();
268            if x.0 >= t {
269                x = x + -t;
270                i += 1 << k;
271            }
272        }
273        i += x.0.round() as i64;
274        if is_neg {
275            i = -i;
276        }
277        i
278    }
279}
280
281impl From<f64> for QuadDouble {
282    fn from(x: f64) -> Self {
283        Self(x, 0., 0., 0.)
284    }
285}
286
287impl Display for QuadDouble {
288    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
289        write!(
290            f,
291            "{}",
292            Decimal::from(self.0)
293                + Decimal::from(self.1)
294                + Decimal::from(self.2)
295                + Decimal::from(self.3)
296        )
297    }
298}
299
300#[derive(Debug, Clone)]
301pub enum ParseDoubleDoubleError {
302    ParseFloatError(ParseFloatError),
303    ParseDecimalError(super::decimal::convert::ParseDecimalError),
304}
305
306impl From<ParseFloatError> for ParseDoubleDoubleError {
307    fn from(e: ParseFloatError) -> Self {
308        Self::ParseFloatError(e)
309    }
310}
311
312impl From<super::decimal::convert::ParseDecimalError> for ParseDoubleDoubleError {
313    fn from(e: super::decimal::convert::ParseDecimalError) -> Self {
314        Self::ParseDecimalError(e)
315    }
316}
317
318impl FromStr for QuadDouble {
319    type Err = ParseDoubleDoubleError;
320    fn from_str(s: &str) -> Result<Self, Self::Err> {
321        let f0: f64 = s.parse()?;
322        let d1 = Decimal::from_str(s)? - Decimal::from(f0);
323        let f1: f64 = d1.to_string().parse()?;
324        let d2 = d1 - Decimal::from(f1);
325        let f2: f64 = d2.to_string().parse()?;
326        let d3 = d2 - Decimal::from(f2);
327        let f3: f64 = d3.to_string().parse()?;
328        Ok(Self::renormalize(f0, f1, f2, f3, 0.))
329    }

Source §

impl QuadDouble

Source

pub fn is_zero(&self) -> bool

Examples found in repository ?

crates/competitive/src/num/quad_double.rs (line 395)

394    pub fn sqrt(self) -> Self {
395        if self.is_zero() {
396            return Self::from(0.);
397        }
398        let x = Self::from(1. / self.0.sqrt());
399        let x = x + x * (Self::from(1.) - self * x * x).div2(2.);
400        let x = x + x * (Self::from(1.) - self * x * x).div2(2.);
401        let x = x + x * (Self::from(1.) - self * x * x).div2(2.);
402        x * self
403    }

Source

pub fn sqrt(self) -> Self

Source

pub fn abs(self) -> Self

Source

fn div2(self, rhs: f64) -> Self

Examples found in repository ?

crates/competitive/src/num/quad_double.rs (line 399)

394    pub fn sqrt(self) -> Self {
395        if self.is_zero() {
396            return Self::from(0.);
397        }
398        let x = Self::from(1. / self.0.sqrt());
399        let x = x + x * (Self::from(1.) - self * x * x).div2(2.);
400        let x = x + x * (Self::from(1.) - self * x * x).div2(2.);
401        let x = x + x * (Self::from(1.) - self * x * x).div2(2.);
402        x * self
403    }