treams_core/special/integrals/
double.rs

1//! Double-double arithmetic: reals `hi + lo` with `|lo|` at most half an ulp of `hi`, and
2//! complex numbers with such parts.
3//!
4//! [`SmallSplitKambe`](super::SmallSplitKambe) recurs the coefficients `c^j / j!` of its
5//! series in this arithmetic and rounds them to double at the end. treams-rs extension.
6
7use crate::Complex;
8
9/// A double-double real `hi + lo`, with `|lo|` at most half an ulp of `hi`.
10#[derive(Clone, Copy, Debug)]
11struct Double {
12    hi: f64,
13    lo: f64,
14}
15
16impl Double {
17    /// `hi + lo` as a double-double.
18    fn normalized(hi: f64, lo: f64) -> Self {
19        let sum = hi + lo;
20        Self {
21            hi: sum,
22            lo: lo - (sum - hi),
23        }
24    }
25
26    /// The exact product of two doubles.
27    fn product(a: f64, b: f64) -> Self {
28        let hi = a * b;
29        Self {
30            hi,
31            lo: a.mul_add(b, -hi),
32        }
33    }
34
35    /// The exact sum of two doubles.
36    fn sum(a: f64, b: f64) -> Self {
37        let hi = a + b;
38        let b_part = hi - a;
39        Self {
40            hi,
41            lo: (a - (hi - b_part)) + (b - b_part),
42        }
43    }
44
45    /// The quotient `self / divisor`, with one correction step of the double quotient.
46    fn divide(self, divisor: Self) -> Self {
47        let quotient = self.hi / divisor.hi;
48        let remainder = self - Self::from(quotient) * divisor;
49        Self::normalized(quotient, remainder.hi / divisor.hi)
50    }
51}
52
53impl From<f64> for Double {
54    fn from(hi: f64) -> Self {
55        Self { hi, lo: 0.0 }
56    }
57}
58
59impl std::ops::Add for Double {
60    type Output = Self;
61    fn add(self, other: Self) -> Self {
62        let sum = Self::sum(self.hi, other.hi);
63        Self::normalized(sum.hi, sum.lo + self.lo + other.lo)
64    }
65}
66
67impl std::ops::Neg for Double {
68    type Output = Self;
69    fn neg(self) -> Self {
70        Self {
71            hi: -self.hi,
72            lo: -self.lo,
73        }
74    }
75}
76
77impl std::ops::Sub for Double {
78    type Output = Self;
79    fn sub(self, other: Self) -> Self {
80        self + -other
81    }
82}
83
84impl std::ops::Mul for Double {
85    type Output = Self;
86    fn mul(self, other: Self) -> Self {
87        let product = Self::product(self.hi, other.hi);
88        let lo = self
89            .hi
90            .mul_add(other.lo, self.lo.mul_add(other.hi, product.lo));
91        Self::normalized(product.hi, lo)
92    }
93}
94
95/// A complex number with double-double parts.
96#[derive(Clone, Copy, Debug)]
97pub(super) struct DoubleComplex {
98    re: Double,
99    im: Double,
100}
101
102impl DoubleComplex {
103    /// `1 / (2 eta^2)` in double-double arithmetic.
104    pub(super) fn half_inverse_square(eta: Complex) -> Self {
105        let (re, im) = (Double::from(eta.re), Double::from(eta.im));
106        // eta^2 = u + i v, and 1 / (2 eta^2) = (u - i v) / (2 (u^2 + v^2)).
107        let (u, v) = (re * re - im * im, Double::from(2.0) * re * im);
108        let scale = Double::from(2.0) * (u * u + v * v);
109        Self {
110            re: u.divide(scale),
111            im: (-v).divide(scale),
112        }
113    }
114
115    /// The quotient of each part by a double.
116    pub(super) fn divide(self, divisor: f64) -> Self {
117        Self {
118            re: self.re.divide(Double::from(divisor)),
119            im: self.im.divide(Double::from(divisor)),
120        }
121    }
122
123    /// The nearest double-precision complex number.
124    pub(super) fn rounded(self) -> Complex {
125        Complex::new(self.re.hi, self.im.hi)
126    }
127}
128
129impl From<Complex> for DoubleComplex {
130    fn from(value: Complex) -> Self {
131        Self {
132            re: Double::from(value.re),
133            im: Double::from(value.im),
134        }
135    }
136}
137
138impl std::ops::Mul for DoubleComplex {
139    type Output = Self;
140    fn mul(self, other: Self) -> Self {
141        Self {
142            re: self.re * other.re - self.im * other.im,
143            im: self.re * other.im + self.im * other.re,
144        }
145    }
146}
treams_core/special/integrals/double.rs

treams_core/special/integrals/
double.rs
