zenith_float_num/
complex_elliptic.rs

1//! Complex elliptic integrals via Carlson \(R_F,R_C,R_D,R_J\).
2//!
3//! Duplication runs in \(\mathbb{C}\) with the principal square root. Parameter
4//! \(m=k^2\); incomplete argument \(x=\sin\varphi\). Not the real series at
5//! \(\lvert z\rvert\).
6
7use crate::common::util::round_p;
8use crate::complex_special::nan_pair;
9use crate::complex_special::pi_c;
10use crate::complex_special::ziv_complex;
11use crate::Consts;
12use crate::Error;
13use crate::ExactComplex;
14use crate::ExactNum;
15use crate::RoundingMode;
16use crate::INF_POS;
17
18/// Same cap as the real Carlson kernels.
19const CARLSON_DUPE_MAX: u32 = 128;
20
21fn rm() -> RoundingMode {
22    RoundingMode::None
23}
24
25fn c_u32(n: u32, p: usize) -> ExactComplex {
26    ExactComplex::from_real(ExactNum::from_u32(n, p), p)
27}
28
29fn is_c_zero(z: &ExactComplex) -> bool {
30    z.re().is_zero() && z.im().is_zero()
31}
32
33fn is_c_one(z: &ExactComplex, p: usize) -> bool {
34    z.im().is_zero() && z.re().cmp(&ExactNum::from_u8(1, p)) == Some(0)
35}
36
37fn c_abs(z: &ExactComplex, p: usize) -> ExactNum {
38    z.abs(p, rm())
39}
40
41fn max_abs(a: &ExactNum, b: &ExactNum, _p: usize) -> ExactNum {
42    if matches!(a.cmp(b), Some(c) if c >= 0) {
43        a.clone()
44    } else {
45        b.clone()
46    }
47}
48
49fn tiny_spread(p: usize) -> ExactNum {
50    ExactNum::from_u8(1, p).ldexp(-((p as i32) / 3 + 16), p, rm())
51}
52
53fn close_enough(dev: &ExactNum, an: &ExactComplex, p: usize) -> bool {
54    let one = ExactNum::from_u8(1, p);
55    let scale = max_abs(&c_abs(an, p), &one, p);
56    let thresh = tiny_spread(p).mul(&scale, p, rm());
57    matches!(dev.cmp(&thresh), Some(c) if c < 0)
58}
59
60fn max_dev3(
61    an: &ExactComplex,
62    x: &ExactComplex,
63    y: &ExactComplex,
64    z: &ExactComplex,
65    p: usize,
66) -> ExactNum {
67    let dx = c_abs(&an.sub(x, p, rm()), p);
68    let dy = c_abs(&an.sub(y, p, rm()), p);
69    let dz = c_abs(&an.sub(z, p, rm()), p);
70    max_abs(&max_abs(&dx, &dy, p), &dz, p)
71}
72
73fn any_nan(args: &[&ExactComplex]) -> bool {
74    args.iter().any(|z| z.is_nan())
75}
76
77/// \(R_C(x,y)=R_F(x,y,y)\).
78fn carlson_rc(x: &ExactComplex, y: &ExactComplex, p: usize, cc: &mut Consts) -> ExactComplex {
79    if any_nan(&[x, y]) {
80        return ExactComplex::new(x.re().clone(), y.im().clone());
81    }
82    if is_c_zero(y) {
83        return nan_pair(Error::InvalidArgument);
84    }
85    let dxy = x.sub(y, p, rm());
86    if close_enough(&c_abs(&dxy, p), x, p) {
87        return c_u32(1, p).div(&x.sqrt(p, rm(), cc), p, rm());
88    }
89    if is_c_zero(x) {
90        let mut hp = pi_c(p, cc);
91        hp = hp.div(&c_u32(2, p), p, rm());
92        return hp.div(&y.sqrt(p, rm(), cc), p, rm());
93    }
94    carlson_rf(x, y, y, p, cc)
95}
96
97/// Symmetric \(R_F(x,y,z)\). At most one argument may be \(0\).
98fn carlson_rf(
99    x0: &ExactComplex,
100    y0: &ExactComplex,
101    z0: &ExactComplex,
102    p: usize,
103    cc: &mut Consts,
104) -> ExactComplex {
105    if any_nan(&[x0, y0, z0]) {
106        return nan_pair(Error::InvalidArgument);
107    }
108    let zeros =
109        usize::from(is_c_zero(x0)) + usize::from(is_c_zero(y0)) + usize::from(is_c_zero(z0));
110    if zeros > 1 {
111        return nan_pair(Error::InvalidArgument);
112    }
113    let four = c_u32(4, p);
114    let three = c_u32(3, p);
115    let mut x = x0.clone();
116    let mut y = y0.clone();
117    let mut z = z0.clone();
118    for _ in 0..CARLSON_DUPE_MAX {
119        let an = x.add(&y, p, rm()).add(&z, p, rm()).div(&three, p, rm());
120        if close_enough(&max_dev3(&an, &x, &y, &z, p), &an, p) {
121            return rf_series(&an, &x, &y, &z, p, cc);
122        }
123        let sx = x.sqrt(p, rm(), cc);
124        let sy = y.sqrt(p, rm(), cc);
125        let sz = z.sqrt(p, rm(), cc);
126        let lam = sx
127            .mul(&sy, p, rm())
128            .add(&sy.mul(&sz, p, rm()), p, rm())
129            .add(&sz.mul(&sx, p, rm()), p, rm());
130        x = x.add(&lam, p, rm()).div(&four, p, rm());
131        y = y.add(&lam, p, rm()).div(&four, p, rm());
132        z = z.add(&lam, p, rm()).div(&four, p, rm());
133    }
134    nan_pair(Error::InvalidArgument)
135}
136
137fn rf_series(
138    an: &ExactComplex,
139    x: &ExactComplex,
140    y: &ExactComplex,
141    z: &ExactComplex,
142    p: usize,
143    cc: &mut Consts,
144) -> ExactComplex {
145    let xx = an.sub(x, p, rm()).div(an, p, rm());
146    let yy = an.sub(y, p, rm()).div(an, p, rm());
147    let zz = an.sub(z, p, rm()).div(an, p, rm());
148    let e2 = xx.mul(&yy, p, rm()).sub(&zz.mul(&zz, p, rm()), p, rm());
149    let e3 = xx.mul(&yy, p, rm()).mul(&zz, p, rm());
150    let e2s = e2.mul(&e2, p, rm());
151    let one = c_u32(1, p);
152    let w = |n: u32| c_u32(n, p);
153    let s = one
154        .sub(&e2.div(&w(10), p, rm()), p, rm())
155        .add(&e3.div(&w(14), p, rm()), p, rm())
156        .add(&e2s.div(&w(24), p, rm()), p, rm())
157        .sub(
158            &w(3)
159                .mul(&e2, p, rm())
160                .mul(&e3, p, rm())
161                .div(&w(44), p, rm()),
162            p,
163            rm(),
164        )
165        .sub(
166            &w(5)
167                .mul(&e2, p, rm())
168                .mul(&e2s, p, rm())
169                .div(&w(208), p, rm()),
170            p,
171            rm(),
172        )
173        .add(
174            &w(3)
175                .mul(&e2s, p, rm())
176                .mul(&e3, p, rm())
177                .div(&w(104), p, rm()),
178            p,
179            rm(),
180        );
181    s.div(&an.sqrt(p, rm(), cc), p, rm())
182}
183
184/// \(R_D(x,y,z)\). \(z\neq 0\); at most one of \(x,y\) may be \(0\).
185fn carlson_rd(
186    x0: &ExactComplex,
187    y0: &ExactComplex,
188    z0: &ExactComplex,
189    p: usize,
190    cc: &mut Consts,
191) -> ExactComplex {
192    if any_nan(&[x0, y0, z0]) {
193        return nan_pair(Error::InvalidArgument);
194    }
195    if is_c_zero(z0) || (is_c_zero(x0) && is_c_zero(y0)) {
196        return nan_pair(Error::InvalidArgument);
197    }
198    let four = c_u32(4, p);
199    let three = c_u32(3, p);
200    let five = c_u32(5, p);
201    let mut x = x0.clone();
202    let mut y = y0.clone();
203    let mut z = z0.clone();
204    let mut sum = ExactComplex::zero(p);
205    let mut fac = c_u32(1, p);
206    for _ in 0..CARLSON_DUPE_MAX {
207        let an = x
208            .add(&y, p, rm())
209            .add(&three.mul(&z, p, rm()), p, rm())
210            .div(&five, p, rm());
211        if close_enough(&max_dev3(&an, &x, &y, &z, p), &an, p) {
212            let series = rd_series(&an, &x, &y, &z, p, cc);
213            return three
214                .mul(&sum, p, rm())
215                .add(&fac.mul(&series, p, rm()), p, rm());
216        }
217        let sx = x.sqrt(p, rm(), cc);
218        let sy = y.sqrt(p, rm(), cc);
219        let sz = z.sqrt(p, rm(), cc);
220        let lam = sx
221            .mul(&sy, p, rm())
222            .add(&sy.mul(&sz, p, rm()), p, rm())
223            .add(&sz.mul(&sx, p, rm()), p, rm());
224        sum = sum.add(
225            &fac.div(&sz.mul(&z.add(&lam, p, rm()), p, rm()), p, rm()),
226            p,
227            rm(),
228        );
229        fac = fac.div(&four, p, rm());
230        x = x.add(&lam, p, rm()).div(&four, p, rm());
231        y = y.add(&lam, p, rm()).div(&four, p, rm());
232        z = z.add(&lam, p, rm()).div(&four, p, rm());
233    }
234    nan_pair(Error::InvalidArgument)
235}
236
237fn rd_series(
238    an: &ExactComplex,
239    x: &ExactComplex,
240    y: &ExactComplex,
241    z: &ExactComplex,
242    p: usize,
243    cc: &mut Consts,
244) -> ExactComplex {
245    let xx = an.sub(x, p, rm()).div(an, p, rm());
246    let yy = an.sub(y, p, rm()).div(an, p, rm());
247    let zz = an.sub(z, p, rm()).div(an, p, rm());
248    let e2 = xx.mul(&yy, p, rm()).sub(&zz.mul(&zz, p, rm()), p, rm());
249    let e3 = xx.mul(&yy, p, rm()).mul(&zz, p, rm());
250    let e2s = e2.mul(&e2, p, rm());
251    let one = c_u32(1, p);
252    let w = |n: u32| c_u32(n, p);
253    let s = one
254        .sub(&w(3).mul(&e2, p, rm()).div(&w(14), p, rm()), p, rm())
255        .add(&e3.div(&w(6), p, rm()), p, rm())
256        .add(&w(9).mul(&e2s, p, rm()).div(&w(88), p, rm()), p, rm())
257        .sub(
258            &w(3)
259                .mul(&e2, p, rm())
260                .mul(&e3, p, rm())
261                .div(&w(22), p, rm()),
262            p,
263            rm(),
264        )
265        .add(
266            &w(9)
267                .mul(&e3, p, rm())
268                .mul(&e3, p, rm())
269                .div(&w(52), p, rm()),
270            p,
271            rm(),
272        )
273        .sub(
274            &w(3)
275                .mul(&e2, p, rm())
276                .mul(&e2s, p, rm())
277                .div(&w(26), p, rm()),
278            p,
279            rm(),
280        );
281    s.div(&an.mul(&an.sqrt(p, rm(), cc), p, rm()), p, rm())
282}
283
284/// \(R_J(x,y,z,p)\). Characteristic \(p\neq 0\); at most one of \(x,y,z\) may be \(0\).
285fn carlson_rj(
286    x0: &ExactComplex,
287    y0: &ExactComplex,
288    z0: &ExactComplex,
289    p0: &ExactComplex,
290    p: usize,
291    cc: &mut Consts,
292) -> ExactComplex {
293    if any_nan(&[x0, y0, z0, p0]) {
294        return nan_pair(Error::InvalidArgument);
295    }
296    if is_c_zero(p0) {
297        return nan_pair(Error::InvalidArgument);
298    }
299    let zeros =
300        usize::from(is_c_zero(x0)) + usize::from(is_c_zero(y0)) + usize::from(is_c_zero(z0));
301    if zeros > 1 {
302        return nan_pair(Error::InvalidArgument);
303    }
304    let two = c_u32(2, p);
305    let three = c_u32(3, p);
306    let four = c_u32(4, p);
307    let five = c_u32(5, p);
308    let mut x = x0.clone();
309    let mut y = y0.clone();
310    let mut z = z0.clone();
311    let mut pv = p0.clone();
312    let mut sum = ExactComplex::zero(p);
313    let mut fac = c_u32(1, p);
314    for _ in 0..CARLSON_DUPE_MAX {
315        let an = x
316            .add(&y, p, rm())
317            .add(&z, p, rm())
318            .add(&two.mul(&pv, p, rm()), p, rm())
319            .div(&five, p, rm());
320        let d4 = max_dev3(&an, &x, &y, &z, p);
321        let dp = c_abs(&an.sub(&pv, p, rm()), p);
322        let dmax = max_abs(&d4, &dp, p);
323        if close_enough(&dmax, &an, p) {
324            let series = rj_series(&an, &x, &y, &z, &pv, p, cc);
325            return three
326                .mul(&sum, p, rm())
327                .add(&fac.mul(&series, p, rm()), p, rm());
328        }
329        let sx = x.sqrt(p, rm(), cc);
330        let sy = y.sqrt(p, rm(), cc);
331        let sz = z.sqrt(p, rm(), cc);
332        let lam = sx
333            .mul(&sy, p, rm())
334            .add(&sy.mul(&sz, p, rm()), p, rm())
335            .add(&sz.mul(&sx, p, rm()), p, rm());
336        let alpha = pv
337            .mul(&sx.add(&sy, p, rm()).add(&sz, p, rm()), p, rm())
338            .add(&sx.mul(&sy, p, rm()).mul(&sz, p, rm()), p, rm());
339        let alpha = alpha.mul(&alpha, p, rm());
340        let pl = pv.add(&lam, p, rm());
341        let beta = pv.mul(&pl, p, rm()).mul(&pl, p, rm());
342        sum = sum.add(
343            &fac.mul(&carlson_rc(&alpha, &beta, p, cc), p, rm()),
344            p,
345            rm(),
346        );
347        fac = fac.div(&four, p, rm());
348        x = x.add(&lam, p, rm()).div(&four, p, rm());
349        y = y.add(&lam, p, rm()).div(&four, p, rm());
350        z = z.add(&lam, p, rm()).div(&four, p, rm());
351        pv = pv.add(&lam, p, rm()).div(&four, p, rm());
352    }
353    nan_pair(Error::InvalidArgument)
354}
355
356fn rj_series(
357    an: &ExactComplex,
358    x: &ExactComplex,
359    y: &ExactComplex,
360    z: &ExactComplex,
361    pv: &ExactComplex,
362    p: usize,
363    cc: &mut Consts,
364) -> ExactComplex {
365    let xx = an.sub(x, p, rm()).div(an, p, rm());
366    let yy = an.sub(y, p, rm()).div(an, p, rm());
367    let zz = an.sub(z, p, rm()).div(an, p, rm());
368    let pp = an.sub(pv, p, rm()).div(an, p, rm());
369    let xyz = xx.mul(&yy, p, rm()).mul(&zz, p, rm());
370    let xy_xz_yz = xx
371        .mul(&yy, p, rm())
372        .add(&xx.mul(&zz, p, rm()), p, rm())
373        .add(&yy.mul(&zz, p, rm()), p, rm());
374    let p2 = pp.mul(&pp, p, rm());
375    let p3 = p2.mul(&pp, p, rm());
376    let two = c_u32(2, p);
377    let three = c_u32(3, p);
378    let e2 = xy_xz_yz.sub(&three.mul(&p2, p, rm()), p, rm());
379    let e3 = xyz
380        .add(&two.mul(&p3, p, rm()), p, rm())
381        .sub(&pp.mul(&xy_xz_yz, p, rm()), p, rm());
382    let e2s = e2.mul(&e2, p, rm());
383    let one = c_u32(1, p);
384    let w = |n: u32| c_u32(n, p);
385    let s = one
386        .sub(&w(3).mul(&e2, p, rm()).div(&w(14), p, rm()), p, rm())
387        .add(&e3.div(&w(6), p, rm()), p, rm())
388        .add(&w(9).mul(&e2s, p, rm()).div(&w(88), p, rm()), p, rm())
389        .sub(
390            &w(3)
391                .mul(&e2, p, rm())
392                .mul(&e3, p, rm())
393                .div(&w(22), p, rm()),
394            p,
395            rm(),
396        )
397        .add(
398            &w(9)
399                .mul(&e3, p, rm())
400                .mul(&e3, p, rm())
401                .div(&w(52), p, rm()),
402            p,
403            rm(),
404        )
405        .sub(
406            &w(3)
407                .mul(&e2, p, rm())
408                .mul(&e2s, p, rm())
409                .div(&w(26), p, rm()),
410            p,
411            rm(),
412        );
413    s.div(&an.mul(&an.sqrt(p, rm(), cc), p, rm()), p, rm())
414}
415
416impl ExactComplex {
417    /// Complete elliptic \(K(m)\), \(m=k^2\). Cut on \([1,+\infty)\). \(m=1\) is \(+\infty\).
418    ///
419    /// # Precision
420    ///
421    /// - Algorithm: Carlson `R_F` in ℂ; `CARLSON_DUPE_MAX = 128`.
422    /// - Bound: identities evaluated outside Ziv (nested Ziv would exhaust `MAX_PREC_RETRY`).
423    /// - MPFR oracle: no.
424    pub fn elliptic_k(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
425        if self.is_nan() {
426            return ExactComplex::new(self.re().clone(), self.im().clone());
427        }
428        let dest = round_p(p);
429        if is_c_one(self, dest) {
430            return ExactComplex::from_real(INF_POS.clone(), dest);
431        }
432        if is_c_zero(self) {
433            let hp = pi_c(dest, cc);
434            return hp.div(&c_u32(2, dest), dest, rm);
435        }
436        ziv_complex(dest, rm, |pw| self.elliptic_k_at(pw, cc))
437    }
438
439    fn elliptic_k_at(&self, p: usize, cc: &mut Consts) -> Self {
440        let zero = ExactComplex::zero(p);
441        let one = ExactComplex::one(p);
442        let om = one.sub(self, p, rm());
443        carlson_rf(&zero, &om, &one, p, cc)
444    }
445
446    /// Complete elliptic \(E(m)\). \(E(1)=1\). Cut of \(K\) inherited through \(1-m\).
447    ///
448    /// # Precision
449    ///
450    /// - Algorithm: Carlson `R_F` / `R_D`; `CARLSON_DUPE_MAX = 128`.
451    /// - Bound: same as [`Self::elliptic_k`].
452    /// - MPFR oracle: no.
453    pub fn elliptic_e_complete(&self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
454        if self.is_nan() {
455            return ExactComplex::new(self.re().clone(), self.im().clone());
456        }
457        let dest = round_p(p);
458        if is_c_zero(self) {
459            let hp = pi_c(dest, cc);
460            return hp.div(&c_u32(2, dest), dest, rm);
461        }
462        if is_c_one(self, dest) {
463            return ExactComplex::one(dest);
464        }
465        ziv_complex(dest, rm, |pw| self.elliptic_e_complete_at(pw, cc))
466    }
467
468    fn elliptic_e_complete_at(&self, p: usize, cc: &mut Consts) -> Self {
469        let zero = ExactComplex::zero(p);
470        let one = ExactComplex::one(p);
471        let three = c_u32(3, p);
472        let om = one.sub(self, p, rm());
473        let rf = carlson_rf(&zero, &om, &one, p, cc);
474        let rd = carlson_rd(&zero, &om, &one, p, cc);
475        rf.sub(&self.div(&three, p, rm()).mul(&rd, p, rm()), p, rm())
476    }
477
478    /// Incomplete \(F(x|m)\), \(x=\sin\varphi\), \(m=k^2\).
479    ///
480    /// # Precision
481    ///
482    /// - Algorithm: Carlson `R_F`; `CARLSON_DUPE_MAX = 128`.
483    /// - MPFR oracle: no.
484    ///
485    /// Carlson \(R_F(1-x^2,1-mx^2,1)\). Cuts when \(1-x^2\) or \(1-mx^2\) lies on
486    /// \((-\infty,0]\) (principal square-root cut). \(F(x,0)=\arcsin x\).
487    pub fn elliptic_f(&self, m: &Self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
488        if self.is_nan() || m.is_nan() {
489            return nan_pair(Error::InvalidArgument);
490        }
491        let dest = round_p(p);
492        if is_c_zero(self) {
493            return ExactComplex::zero(dest);
494        }
495        if is_c_zero(m) {
496            return self.asin(dest, rm, cc);
497        }
498        if is_c_one(m, dest) {
499            return self.atanh(dest, rm, cc);
500        }
501        ziv_complex(dest, rm, |pw| self.elliptic_f_at(m, pw, cc))
502    }
503
504    fn elliptic_f_at(&self, m: &Self, p: usize, cc: &mut Consts) -> Self {
505        let one = ExactComplex::one(p);
506        let x2 = self.mul(self, p, rm());
507        let a = one.sub(&x2, p, rm());
508        let b = one.sub(&m.mul(&x2, p, rm()), p, rm());
509        let rf = carlson_rf(&a, &b, &one, p, cc);
510        self.mul(&rf, p, rm())
511    }
512
513    /// Incomplete \(E(x|m)\). Same \(x,m\) convention as [`Self::elliptic_f`].
514    ///
515    /// # Precision
516    ///
517    /// - Algorithm: Carlson `R_F` / `R_D`; `CARLSON_DUPE_MAX = 128`.
518    /// - MPFR oracle: no.
519    /// Cuts as for \(F\). \(E(x,0)=\arcsin x\); \(E(x,1)=x\).
520    pub fn elliptic_e(&self, m: &Self, p: usize, rm: RoundingMode, cc: &mut Consts) -> Self {
521        if self.is_nan() || m.is_nan() {
522            return nan_pair(Error::InvalidArgument);
523        }
524        let dest = round_p(p);
525        if is_c_zero(self) {
526            return ExactComplex::zero(dest);
527        }
528        if is_c_zero(m) {
529            return self.asin(dest, rm, cc);
530        }
531        if is_c_one(m, dest) {
532            return self.clone();
533        }
534        ziv_complex(dest, rm, |pw| self.elliptic_e_at(m, pw, cc))
535    }
536
537    fn elliptic_e_at(&self, m: &Self, p: usize, cc: &mut Consts) -> Self {
538        let one = ExactComplex::one(p);
539        let three = c_u32(3, p);
540        let x2 = self.mul(self, p, rm());
541        let a = one.sub(&x2, p, rm());
542        let b = one.sub(&m.mul(&x2, p, rm()), p, rm());
543        let rf = carlson_rf(&a, &b, &one, p, cc);
544        let rd = carlson_rd(&a, &b, &one, p, cc);
545        self.mul(&rf, p, rm()).sub(
546            &m.mul(self, p, rm())
547                .mul(&x2, p, rm())
548                .div(&three, p, rm())
549                .mul(&rd, p, rm()),
550            p,
551            rm(),
552        )
553    }
554
555    /// Complete \(\Pi(n,m)\). `self` is \(n\). \(\Pi(0,m)=K(m)\). Pole at \(n=1\).
556    ///
557    /// # Precision
558    ///
559    /// - Algorithm: Carlson `R_J`; `CARLSON_DUPE_MAX = 128`.
560    /// - MPFR oracle: no.
561    pub fn elliptic_pi_complete(
562        &self,
563        m: &Self,
564        p: usize,
565        rm: RoundingMode,
566        cc: &mut Consts,
567    ) -> Self {
568        if self.is_nan() || m.is_nan() {
569            return nan_pair(Error::InvalidArgument);
570        }
571        let dest = round_p(p);
572        if is_c_one(self, dest) {
573            return nan_pair(Error::InvalidArgument);
574        }
575        if is_c_zero(self) {
576            return m.elliptic_k(dest, rm, cc);
577        }
578        ziv_complex(dest, rm, |pw| self.elliptic_pi_complete_at(m, pw, cc))
579    }
580
581    fn elliptic_pi_complete_at(&self, m: &Self, p: usize, cc: &mut Consts) -> Self {
582        let zero = ExactComplex::zero(p);
583        let one = ExactComplex::one(p);
584        let three = c_u32(3, p);
585        let om = one.sub(m, p, rm());
586        let on = one.sub(self, p, rm());
587        let rf = carlson_rf(&zero, &om, &one, p, cc);
588        let rj = carlson_rj(&zero, &om, &one, &on, p, cc);
589        rf.add(&self.div(&three, p, rm()).mul(&rj, p, rm()), p, rm())
590    }
591
592    /// Incomplete \(\Pi(n;x|m)\). `self` is \(n\).
593    ///
594    /// # Precision
595    ///
596    /// - Algorithm: Carlson `R_J`; `CARLSON_DUPE_MAX = 128`.
597    /// - MPFR oracle: no.
598    ///
599    /// Cuts when \(1-x^2\), \(1-mx^2\), or \(1-nx^2\) meets the Carlson cut
600    /// \((-\infty,0]\). \(\Pi(0;x|m)=F(x|m)\).
601    pub fn elliptic_pi(
602        &self,
603        x: &Self,
604        m: &Self,
605        p: usize,
606        rm: RoundingMode,
607        cc: &mut Consts,
608    ) -> Self {
609        if self.is_nan() || x.is_nan() || m.is_nan() {
610            return nan_pair(Error::InvalidArgument);
611        }
612        let dest = round_p(p);
613        if is_c_zero(x) {
614            return ExactComplex::zero(dest);
615        }
616        if is_c_zero(self) {
617            return x.elliptic_f(m, dest, rm, cc);
618        }
619        ziv_complex(dest, rm, |pw| self.elliptic_pi_at(x, m, pw, cc))
620    }
621
622    fn elliptic_pi_at(&self, x: &Self, m: &Self, p: usize, cc: &mut Consts) -> Self {
623        let one = ExactComplex::one(p);
624        let three = c_u32(3, p);
625        let x2 = x.mul(x, p, rm());
626        let a = one.sub(&x2, p, rm());
627        let b = one.sub(&m.mul(&x2, p, rm()), p, rm());
628        let pv = one.sub(&self.mul(&x2, p, rm()), p, rm());
629        let rf = carlson_rf(&a, &b, &one, p, cc);
630        let rj = carlson_rj(&a, &b, &one, &pv, p, cc);
631        x.mul(&rf, p, rm()).add(
632            &self
633                .mul(x, p, rm())
634                .mul(&x2, p, rm())
635                .div(&three, p, rm())
636                .mul(&rj, p, rm()),
637            p,
638            rm(),
639        )
640    }
641}
642
643#[cfg(test)]
644mod tests {
645    use super::*;
646
647    fn near(a: &ExactNum, b: &ExactNum, p: usize, slack: i32) -> bool {
648        let d = a.sub(b, p, RoundingMode::None).abs();
649        d.is_zero() || d.exponent().unwrap_or(0) < -((p as i32) - slack)
650    }
651
652    fn tiny(x: &ExactNum, p: usize) -> bool {
653        x.is_zero() || x.exponent().unwrap_or(0) < -((p as i32) / 4)
654    }
655
656    #[test]
657    fn complex_elliptic_golds() {
658        let p = 256;
659        let r = RoundingMode::ToEven;
660        let mut cc = Consts::new().unwrap();
661        let z0 = ExactComplex::zero(p);
662        let half = ExactNum::from_u8(1, p).div(&ExactNum::from_u8(2, p), p, r);
663        let mh = ExactComplex::from_real(half.clone(), p);
664
665        let k0 = z0.elliptic_k(p, r, &mut cc);
666        let mut hp = cc.pi(p, r);
667        hp = hp.div(&ExactNum::from_u8(2, p), p, r);
668        assert!(near(k0.re(), &hp, p, 40), "K(0) re");
669        assert!(tiny(k0.im(), p), "K(0) im");
670
671        let e0 = z0.elliptic_e_complete(p, r, &mut cc);
672        assert!(near(e0.re(), &hp, p, 40), "E(0) re");
673        assert!(tiny(e0.im(), p), "E(0) im");
674
675        let k_real = half.elliptic_k(p, r, &mut cc);
676        let k_c = mh.elliptic_k(p, r, &mut cc);
677        assert!(near(k_c.re(), &k_real, p, 40), "K(1/2) matches real");
678        assert!(tiny(k_c.im(), p), "K(1/2) im");
679
680        let e1 = ExactComplex::one(p).elliptic_e_complete(p, r, &mut cc);
681        assert!(near(e1.re(), &ExactNum::from_u8(1, p), p, 40), "E(1)");
682        assert!(tiny(e1.im(), p));
683
684        let k1 = ExactComplex::one(p).elliptic_k(p, r, &mut cc);
685        assert!(k1.re().is_inf_pos(), "K(1)=+∞");
686
687        let x = mh.clone();
688        let f0 = x.elliptic_f(&z0, p, r, &mut cc);
689        let asin = x.asin(p, r, &mut cc);
690        assert!(near(f0.re(), asin.re(), p, 8), "F(x,0)=arcsin");
691        assert!(tiny(f0.im(), p) && tiny(asin.im(), p));
692
693        let n0 = z0.elliptic_pi_complete(&mh, p, r, &mut cc);
694        assert!(near(n0.re(), k_c.re(), p, 40), "Π(0,m)=K(m)");
695
696        // Legendre: E(m)K(1-m)+E(1-m)K(m)-K(m)K(1-m)=π/2
697        let m = ExactComplex::new(
698            ExactNum::from_u8(3, p).div(&ExactNum::from_u8(10, p), p, r),
699            ExactNum::from_u8(1, p).div(&ExactNum::from_u8(10, p), p, r),
700        );
701        let mp = ExactComplex::one(p).sub(&m, p, r);
702        let km = m.elliptic_k(p, r, &mut cc);
703        let kp = mp.elliptic_k(p, r, &mut cc);
704        let em = m.elliptic_e_complete(p, r, &mut cc);
705        let ep = mp.elliptic_e_complete(p, r, &mut cc);
706        let lhs = em
707            .mul(&kp, p, r)
708            .add(&ep.mul(&km, p, r), p, r)
709            .sub(&km.mul(&kp, p, r), p, r);
710        let want = ExactComplex::from_real(hp, p);
711        assert!(
712            near(lhs.re(), want.re(), p, 30),
713            "Legendre re {:?}",
714            lhs.re()
715        );
716        assert!(
717            tiny(lhs.im(), p) || near(lhs.im(), want.im(), p, 30),
718            "Legendre im"
719        );
720
721        // Cut of K on [1,+∞): K(2+εi) and K(2-εi) are conjugates, not equal.
722        let two = ExactNum::from_u8(2, p);
723        let eps = ExactNum::from_u8(1, p).ldexp(-40, p, RoundingMode::None);
724        let above = ExactComplex::new(two.clone(), eps.clone());
725        let below = ExactComplex::new(two, eps.neg());
726        let ka = above.elliptic_k(p, r, &mut cc);
727        let kb = below.elliptic_k(p, r, &mut cc);
728        assert!(near(ka.re(), kb.re(), p, 20), "K cut Re");
729        assert!(ka.im().is_positive() != kb.im().is_positive() || !tiny(ka.im(), p));
730        assert!(near(&ka.im().abs(), &kb.im().abs(), p, 20), "K cut |Im|");
731        assert_ne!(ka.im().cmp(kb.im()), Some(0));
732
733        // Incomplete F cut: x=3/2, m=1/2, ±ε imag on x.
734        let three_half = ExactNum::from_u8(3, p).div(&ExactNum::from_u8(2, p), p, r);
735        let xa = ExactComplex::new(three_half.clone(), eps.clone());
736        let xb = ExactComplex::new(three_half, eps.neg());
737        let fa = xa.elliptic_f(&mh, p, r, &mut cc);
738        let fb = xb.elliptic_f(&mh, p, r, &mut cc);
739        assert!(!fa.is_nan() && !fb.is_nan(), "F cut defined");
740        assert_ne!(fa.im().cmp(fb.im()), Some(0), "F cut Im differs");
741    }
742}
zenith_float_num/complex_elliptic.rs

zenith_float_num/
complex_elliptic.rs
