zenith_float_num/
complex_ei.rs

1//! Complex `Ei`, `Si`, `Ci`, `li`, and Fresnel integrals.
2
3use crate::common::util::round_p;
4use crate::complex_special::half_c;
5use crate::complex_special::nan_pair;
6use crate::complex_special::neg_c;
7use crate::complex_special::pi_c;
8use crate::complex_special::series_term_cap;
9use crate::complex_special::term_negligible;
10use crate::complex_special::two_c;
11use crate::complex_special::ziv_complex;
12use crate::Consts;
13use crate::Error;
14use crate::ExactComplex;
15use crate::ExactNum;
16use crate::RoundingMode;
17
18/// Use the power series for \(\mathrm{Ei}\) when \(\lvert z\rvert\) is below this.
19const EI_SERIES_THRESHOLD: u32 = 16;
20
21fn abs_below(z: &ExactComplex, bound: u32, p: usize) -> bool {
22    let a = z.abs(p, RoundingMode::None);
23    let b = ExactNum::from_u32(bound, p);
24    matches!(a.cmp(&b), Some(c) if c < 0)
25}
26
27fn use_ei_series(z: &ExactComplex, dest_p: usize) -> bool {
28    if abs_below(z, EI_SERIES_THRESHOLD, dest_p) {
29        return true;
30    }
31    let az = z.abs(dest_p, RoundingMode::None);
32    let az2 = az.mul(&az, dest_p, RoundingMode::None);
33    let thresh = ExactNum::from_u32(dest_p.min(u32::MAX as usize) as u32, dest_p);
34    matches!(az2.cmp(&thresh), Some(c) if c < 0)
35}
36
37impl ExactComplex {
38    /// Exponential integral \(\mathrm{Ei}(z)\). Cut on \((-\infty,0]\); pole at \(0\) → NaN.
39    ///
40    /// # Precision
41    ///
42    /// - Algorithm: power series for `|z| < EI_SERIES_THRESHOLD` (`16`); asymptotic otherwise.
43    /// - Bound: Ziv on each part (`MAX_PREC_RETRY`).
44    /// - MPFR oracle: no.
45    pub fn ei(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
46        if self.is_nan() {
47            return ExactComplex::new(self.re().clone(), self.im().clone());
48        }
49        if self.re().is_zero() && self.im().is_zero() {
50            return nan_pair(Error::InvalidArgument);
51        }
52        let dest = round_p(p);
53        ziv_complex(dest, rm, |pw| self.ei_at(pw, dest, cc))
54    }
55
56    /// Sine integral \(\mathrm{Si}(z)=(E_i(iz)-E_i(-iz))/(2i)-\pi/2\).
57    ///
58    /// # Precision
59    ///
60    /// - Algorithm: via [`Self::ei`]; inherits `EI_SERIES_THRESHOLD = 16`.
61    /// - Bound: Ziv on each part (`MAX_PREC_RETRY`).
62    /// - MPFR oracle: no.
63    pub fn si(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
64        if self.is_nan() {
65            return ExactComplex::new(self.re().clone(), self.im().clone());
66        }
67        let dest = round_p(p);
68        ziv_complex(dest, rm, |pw| self.si_at(pw, dest, cc))
69    }
70
71    /// Cosine integral \(\mathrm{Ci}(z)\). Pole at \(0\) → NaN. Inherits the \(\mathrm{Ei}\) cut.
72    ///
73    /// # Precision
74    ///
75    /// - Algorithm: via [`Self::ei`]; `EI_SERIES_THRESHOLD = 16`.
76    /// - Bound: Ziv on each part (`MAX_PREC_RETRY`).
77    /// - MPFR oracle: no.
78    pub fn ci(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
79        if self.is_nan() {
80            return ExactComplex::new(self.re().clone(), self.im().clone());
81        }
82        if self.re().is_zero() && self.im().is_zero() {
83            return nan_pair(Error::InvalidArgument);
84        }
85        let dest = round_p(p);
86        ziv_complex(dest, rm, |pw| self.ci_at(pw, dest, cc))
87    }
88
89    /// Logarithmic integral \(\mathrm{li}(z)=\mathrm{Ei}(\ln z)\). Cut on \((-\infty,1]\); pole at \(1\) → NaN.
90    ///
91    /// # Precision
92    ///
93    /// - Algorithm: [`Self::ei`] of `ln z`.
94    /// - Bound: Ziv on each part (`MAX_PREC_RETRY`).
95    /// - MPFR oracle: no.
96    pub fn li(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
97        if self.is_nan() {
98            return ExactComplex::new(self.re().clone(), self.im().clone());
99        }
100        if self.re().is_zero() && self.im().is_zero() {
101            return nan_pair(Error::InvalidArgument);
102        }
103        let dest = round_p(p);
104        ziv_complex(dest, rm, |pw| self.li_at(pw, dest, cc))
105    }
106
107    /// Fresnel sine integral \(S(z)\). Entire.
108    ///
109    /// # Precision
110    ///
111    /// - Algorithm: via complex `erf`; Ziv on each part (`MAX_PREC_RETRY`).
112    /// - MPFR oracle: no.
113    pub fn fresnel_s(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
114        if self.is_nan() {
115            return ExactComplex::new(self.re().clone(), self.im().clone());
116        }
117        let dest = round_p(p);
118        ziv_complex(dest, rm, |pw| self.fresnel_s_at(pw, cc))
119    }
120
121    /// Fresnel cosine integral \(C(z)\). Entire.
122    ///
123    /// # Precision
124    ///
125    /// - Algorithm: via complex `erf`; Ziv on each part (`MAX_PREC_RETRY`).
126    /// - MPFR oracle: no.
127    pub fn fresnel_c(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
128        if self.is_nan() {
129            return ExactComplex::new(self.re().clone(), self.im().clone());
130        }
131        let dest = round_p(p);
132        ziv_complex(dest, rm, |pw| self.fresnel_c_at(pw, cc))
133    }
134
135    fn ei_at(&self, work_p: usize, dest_p: usize, cc: &mut Consts) -> Self {
136        if use_ei_series(self, dest_p) {
137            self.ei_series(work_p, cc)
138        } else {
139            self.ei_asymptotic(work_p, cc)
140        }
141    }
142
143    /// \(\mathrm{Ei}(z)=\gamma+\mathrm{Ln}\,z+\sum z^n/(n\cdot n!)\). Principal \(\mathrm{Ln}\).
144    fn ei_series(&self, p: usize, cc: &mut Consts) -> Self {
145        let g = ExactComplex::from_real(cc.euler_gamma(p, RoundingMode::None), p);
146        let lnz = self.ln(p, RoundingMode::None, cc);
147        let mut term = self.clone();
148        let mut sum = term.clone();
149        for n in 2..=series_term_cap(p) {
150            let nw = ExactComplex::from_real(ExactNum::from_u32(n as u32, p), p);
151            term = term
152                .mul(self, p, RoundingMode::None)
153                .div(&nw, p, RoundingMode::None);
154            let piece = term.div(&nw, p, RoundingMode::None);
155            sum = sum.add(&piece, p, RoundingMode::None);
156            if term_negligible(&piece, p) {
157                break;
158            }
159        }
160        g.add(&lnz, p, RoundingMode::None)
161            .add(&sum, p, RoundingMode::None)
162    }
163
164    /// \(\mathrm{Ei}(z)\sim e^z/z\sum k!/z^k\), stopped at the smallest term.
165    fn ei_asymptotic(&self, p: usize, cc: &mut Consts) -> Self {
166        let pre = self
167            .exp(p, RoundingMode::None, cc)
168            .div(self, p, RoundingMode::None);
169        let mut term = ExactComplex::one(p);
170        let mut s = term.clone();
171        let mut prev_e = i32::MIN;
172        for k in 1..=series_term_cap(p) {
173            let kk = ExactComplex::from_real(ExactNum::from_u32(k as u32, p), p);
174            term = term
175                .mul(&kk, p, RoundingMode::None)
176                .div(self, p, RoundingMode::None);
177            let e = term
178                .abs(p, RoundingMode::None)
179                .exponent()
180                .unwrap_or(i32::MIN);
181            if k > 1 && e > prev_e {
182                break;
183            }
184            prev_e = e;
185            s = s.add(&term, p, RoundingMode::None);
186            if term_negligible(&term, p) {
187                break;
188            }
189        }
190        pre.mul(&s, p, RoundingMode::None)
191    }
192
193    fn si_at(&self, work_p: usize, dest_p: usize, cc: &mut Consts) -> Self {
194        if self.re().is_zero() && self.im().is_zero() {
195            return ExactComplex::zero(work_p);
196        }
197        // \(\mathrm{Si}\) is odd; the \(\mathrm{Ei}\) combination minus \(\pi/2\) is the
198        // principal value for \(\mathrm{Re}\,z\ge 0\).
199        if self.re().is_negative() || (self.re().is_zero() && self.im().is_negative()) {
200            return neg_c(&neg_c(self).si_at(work_p, dest_p, cc));
201        }
202        let iz = ExactComplex::i(work_p).mul(self, work_p, RoundingMode::None);
203        let e_plus = iz.ei_at(work_p, dest_p, cc);
204        let e_minus = neg_c(&iz).ei_at(work_p, dest_p, cc);
205        let two_i = two_c(work_p).mul(&ExactComplex::i(work_p), work_p, RoundingMode::None);
206        let half_pi = pi_c(work_p, cc).mul(&half_c(work_p), work_p, RoundingMode::None);
207        e_plus
208            .sub(&e_minus, work_p, RoundingMode::None)
209            .div(&two_i, work_p, RoundingMode::None)
210            .sub(&half_pi, work_p, RoundingMode::None)
211    }
212
213    fn ci_at(&self, work_p: usize, dest_p: usize, cc: &mut Consts) -> Self {
214        let iz = ExactComplex::i(work_p).mul(self, work_p, RoundingMode::None);
215        let e_plus = iz.ei_at(work_p, dest_p, cc);
216        let e_minus = neg_c(&iz).ei_at(work_p, dest_p, cc);
217        neg_c(&e_plus.add(&e_minus, work_p, RoundingMode::None).div(
218            &two_c(work_p),
219            work_p,
220            RoundingMode::None,
221        ))
222    }
223
224    fn li_at(&self, work_p: usize, dest_p: usize, cc: &mut Consts) -> Self {
225        let one = ExactComplex::one(work_p);
226        if self.im().is_zero() && self.re().cmp(one.re()) == Some(0) {
227            return nan_pair(Error::InvalidArgument);
228        }
229        self.ln(work_p, RoundingMode::None, cc)
230            .ei_at(work_p, dest_p, cc)
231    }
232
233    fn fresnel_pair(&self, p: usize, cc: &mut Consts) -> (Self, Self) {
234        let sqrt_pi =
235            ExactComplex::from_real(cc.pi(p, RoundingMode::None).sqrt(p, RoundingMode::None), p);
236        let scale = sqrt_pi
237            .mul(self, p, RoundingMode::None)
238            .mul(&half_c(p), p, RoundingMode::None);
239        let one = ExactComplex::one(p);
240        let i = ExactComplex::i(p);
241        let one_p_i = one.add(&i, p, RoundingMode::None);
242        let one_m_i = one.sub(&i, p, RoundingMode::None);
243        let erf_m = one_m_i
244            .mul(&scale, p, RoundingMode::None)
245            .erf(p, RoundingMode::None, cc);
246        let erf_p = one_p_i
247            .mul(&scale, p, RoundingMode::None)
248            .erf(p, RoundingMode::None, cc);
249        let c_plus_is =
250            one_p_i
251                .mul(&half_c(p), p, RoundingMode::None)
252                .mul(&erf_m, p, RoundingMode::None);
253        let c_minus_is =
254            one_m_i
255                .mul(&half_c(p), p, RoundingMode::None)
256                .mul(&erf_p, p, RoundingMode::None);
257        let c = c_plus_is.add(&c_minus_is, p, RoundingMode::None).mul(
258            &half_c(p),
259            p,
260            RoundingMode::None,
261        );
262        let two_i = two_c(p).mul(&i, p, RoundingMode::None);
263        let s =
264            c_plus_is
265                .sub(&c_minus_is, p, RoundingMode::None)
266                .div(&two_i, p, RoundingMode::None);
267        (s, c)
268    }
269
270    fn fresnel_s_at(&self, p: usize, cc: &mut Consts) -> Self {
271        self.fresnel_pair(p, cc).0
272    }
273
274    fn fresnel_c_at(&self, p: usize, cc: &mut Consts) -> Self {
275        self.fresnel_pair(p, cc).1
276    }
277}
278
279#[cfg(test)]
280mod tests {
281    use super::*;
282    use crate::complex_special::neg_c;
283    use crate::complex_special::pi_c;
284    use crate::complex_special::two_c;
285
286    fn near(a: &ExactNum, b: &ExactNum, p: usize) -> bool {
287        let d = a.sub(b, p, RoundingMode::None).abs();
288        d.is_zero() || d.exponent().is_some_and(|e| e < -((p as i32) / 4))
289    }
290
291    fn cnear(a: &ExactComplex, b: &ExactComplex, p: usize) -> bool {
292        near(a.re(), b.re(), p) && near(a.im(), b.im(), p)
293    }
294
295    fn cnear_bits(a: &ExactComplex, b: &ExactComplex, p: usize, slack: i32) -> bool {
296        let dr = a.re().sub(b.re(), p, RoundingMode::None).abs();
297        let di = a.im().sub(b.im(), p, RoundingMode::None).abs();
298        (dr.is_zero() || dr.exponent().is_some_and(|e| e < -((p as i32) / slack)))
299            && (di.is_zero() || di.exponent().is_some_and(|e| e < -((p as i32) / slack)))
300    }
301
302    fn tiny(x: &ExactNum, p: usize) -> bool {
303        x.is_zero() || x.exponent().is_some_and(|e| e < -((p as i32) / 4))
304    }
305
306    #[test]
307    fn test_complex_ei_golds() {
308        let p = 256;
309        let rm = RoundingMode::ToEven;
310        let mut cc = Consts::new().unwrap();
311
312        let one = ExactComplex::one(p);
313        let ei1 = one.ei(p, rm, &mut cc);
314        let r1 = ExactNum::from_u8(1, p).ei(p, rm, &mut cc);
315        assert!(near(ei1.re(), &r1, p));
316        assert!(tiny(ei1.im(), p));
317
318        let z0 = ExactComplex::zero(p);
319        let s0 = z0.si(p, rm, &mut cc);
320        assert!(tiny(s0.re(), p) && tiny(s0.im(), p));
321        assert!(z0.ci(p, rm, &mut cc).is_nan());
322
323        let z = ExactComplex::new(ExactNum::from_u8(1, p), half_c(p).re().clone());
324        let sz = z.si(p, rm, &mut cc);
325        assert!(cnear(&sz, &neg_c(&neg_c(&z).si(p, rm, &mut cc)), p));
326
327        let fs0 = z0.fresnel_s(p, rm, &mut cc);
328        let fc0 = z0.fresnel_c(p, rm, &mut cc);
329        assert!(tiny(fs0.re(), p) && tiny(fs0.im(), p));
330        assert!(tiny(fc0.re(), p) && tiny(fc0.im(), p));
331
332        let h = ExactNum::from_u8(2, p).powsi(-((p as isize) / 8), p, rm);
333        let hc = ExactComplex::from_real(h, p);
334        let num_ei =
335            z.add(&hc, p, rm)
336                .ei(p, rm, &mut cc)
337                .sub(&z.sub(&hc, p, rm).ei(p, rm, &mut cc), p, rm);
338        let deriv_ei = num_ei.div(&hc.mul(&two_c(p), p, rm), p, rm);
339        let expect_ei = z.exp(p, rm, &mut cc).div(&z, p, rm);
340        assert!(cnear_bits(&deriv_ei, &expect_ei, p, 8));
341
342        let num_si =
343            z.add(&hc, p, rm)
344                .si(p, rm, &mut cc)
345                .sub(&z.sub(&hc, p, rm).si(p, rm, &mut cc), p, rm);
346        let deriv_si = num_si.div(&hc.mul(&two_c(p), p, rm), p, rm);
347        let expect_si = z.sin(p, rm, &mut cc).div(&z, p, rm);
348        assert!(cnear_bits(&deriv_si, &expect_si, p, 8));
349
350        let e = ExactComplex::from_real(cc.e(p, rm), p);
351        let lie = e.li(p, rm, &mut cc);
352        let rli = cc.e(p, rm).li(p, rm, &mut cc);
353        assert!(near(lie.re(), &rli, p));
354        assert!(tiny(lie.im(), p));
355
356        let above = ExactComplex::new(ExactNum::from_i8(-1, p), ExactNum::new(p));
357        let below = ExactComplex::new(ExactNum::from_i8(-1, p), ExactNum::new(p).neg());
358        assert!(above.im().is_positive());
359        assert!(below.im().is_negative());
360        let jump = above
361            .ei(p, rm, &mut cc)
362            .sub(&below.ei(p, rm, &mut cc), p, rm);
363        let two_pi_i = two_c(p)
364            .mul(&pi_c(p, &mut cc), p, rm)
365            .mul(&ExactComplex::i(p), p, rm);
366        assert!(cnear_bits(&jump, &two_pi_i, p, 8));
367
368        let nan = ExactComplex::new(crate::NAN.clone(), ExactNum::new(p));
369        assert!(nan.ei(p, rm, &mut cc).is_nan());
370    }
371}
zenith_float_num/complex_ei.rs

zenith_float_num/
complex_ei.rs
