zenith_float_num/
complex_hypergeom.rs

1//! Complex Gaussian \({}_2F_1(a,b;c;z)\).
2//!
3//! Series, Euler, Pfaff, and Kummer / linear transformations in \(\mathbb{C}\).
4//! Not the real series at \(\lvert z\rvert\).
5
6use crate::common::util::round_p;
7use crate::complex_special::is_nonpos_integer;
8use crate::complex_special::nan_pair;
9use crate::complex_special::neg_c;
10use crate::complex_special::term_negligible;
11use crate::complex_special::ziv_complex;
12use crate::Consts;
13use crate::Error;
14use crate::ExactComplex;
15use crate::ExactNum;
16use crate::RoundingMode;
17use crate::WORD_BIT_SIZE;
18
19/// Series term cap (named; same order as the real kernel).
20const HYPERGEOM_SERIES_MAX_TERMS: u32 = 10_000;
21
22/// Linear-transform attempts before `NaN`.
23const HYPERGEOM_TRANSFORM_MAX: u32 = 8;
24
25fn rm() -> RoundingMode {
26    RoundingMode::None
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_eq(a: &ExactComplex, b: &ExactComplex) -> bool {
38    a.re().cmp(b.re()) == Some(0) && a.im().cmp(b.im()) == Some(0)
39}
40
41fn abs_lt_one(z: &ExactComplex, p: usize) -> bool {
42    let az = z.abs(p, rm());
43    matches!(az.cmp(&ExactNum::from_u8(1, p)), Some(c) if c < 0)
44}
45
46fn re_lt_half(z: &ExactComplex, p: usize) -> bool {
47    let half = ExactNum::from_u8(1, p).div(&ExactNum::from_u8(2, p), p, rm());
48    matches!(z.re().cmp(&half), Some(c) if c < 0)
49}
50
51fn re_positive_strict(z: &ExactComplex) -> bool {
52    z.re().is_positive() && !z.re().is_zero()
53}
54
55fn any_nan(args: &[&ExactComplex]) -> bool {
56    args.iter().any(|z| z.is_nan())
57}
58
59fn as_i32_real_int(z: &ExactComplex, p: usize) -> Option<i32> {
60    if !z.im().is_zero() || !z.re().is_int() {
61        return None;
62    }
63    if z.re().is_zero() {
64        return Some(0);
65    }
66    for n in 1i32..=10_000 {
67        let w = ExactNum::from_u32(n as u32, p);
68        if z.re().abs().cmp(&w) == Some(0) {
69            return Some(if z.re().is_negative() { -n } else { n });
70        }
71    }
72    None
73}
74
75fn terminating_neg_int(v: &ExactComplex, p: usize) -> bool {
76    as_i32_real_int(v, p).is_some_and(|n| n <= 0)
77}
78
79fn terminating_after(v: &ExactComplex, next_n: u32, p: usize) -> bool {
80    as_i32_real_int(v, p).is_some_and(|m| m <= 0 && next_n as i32 > -m)
81}
82
83fn pole_c(c: &ExactComplex, a: &ExactComplex, b: &ExactComplex, p: usize) -> bool {
84    let Some(cn) = as_i32_real_int(c, p) else {
85        return false;
86    };
87    if cn > 0 {
88        return false;
89    }
90    let pole_k = -cn;
91    let stop_a = as_i32_real_int(a, p).filter(|&n| n <= 0);
92    let stop_b = as_i32_real_int(b, p).filter(|&n| n <= 0);
93    match (stop_a, stop_b) {
94        (Some(sa), _) if -sa < pole_k => false,
95        (_, Some(sb)) if -sb < pole_k => false,
96        _ => true,
97    }
98}
99
100fn gamma_fixed(z: &ExactComplex, p: usize, cc: &mut Consts) -> ExactComplex {
101    if is_nonpos_integer(z) {
102        return nan_pair(Error::InvalidArgument);
103    }
104    z.gamma_at(p, cc)
105}
106
107fn hypergeom_series(
108    a: &ExactComplex,
109    b: &ExactComplex,
110    c: &ExactComplex,
111    z: &ExactComplex,
112    p: usize,
113) -> ExactComplex {
114    let one = ExactComplex::one(p);
115    let mut term = ExactComplex::one(p);
116    let mut sum = ExactComplex::one(p);
117    let tiny = ExactNum::from_u8(1, p).ldexp(-((p as i32) - 8), p, rm());
118    let n_max = (p.saturating_add(WORD_BIT_SIZE).saturating_add(32))
119        .min(HYPERGEOM_SERIES_MAX_TERMS as usize);
120    for n in 0..n_max {
121        if n > 0 && term_negligible(&term, p) {
122            break;
123        }
124        let nn = ExactComplex::from_real(ExactNum::from_u32(n as u32, p), p);
125        let den = c.add(&nn, p, rm()).mul(&one.add(&nn, p, rm()), p, rm());
126        if matches!(den.abs(p, rm()).cmp(&tiny), Some(c) if c < 0) {
127            return nan_pair(Error::InvalidArgument);
128        }
129        term = term
130            .mul(&a.add(&nn, p, rm()), p, rm())
131            .mul(&b.add(&nn, p, rm()), p, rm())
132            .div(&den, p, rm())
133            .mul(z, p, rm());
134        sum = sum.add(&term, p, rm());
135        if terminating_after(a, (n + 1) as u32, p) || terminating_after(b, (n + 1) as u32, p) {
136            sum.set_inexact(false);
137            return sum;
138        }
139    }
140    sum
141}
142
143/// Kummer: \({}_2F_1(a,b;c;1)=\Gamma(c)\Gamma(c-a-b)/(\Gamma(c-a)\Gamma(c-b))\)
144/// when \(\mathrm{Re}(c-a-b)>0\).
145fn kummer_z_one(
146    a: &ExactComplex,
147    b: &ExactComplex,
148    c: &ExactComplex,
149    p: usize,
150    cc: &mut Consts,
151) -> ExactComplex {
152    let cab = c.sub(a, p, rm()).sub(b, p, rm());
153    if !re_positive_strict(&cab) {
154        return nan_pair(Error::InvalidArgument);
155    }
156    let gc = gamma_fixed(c, p, cc);
157    let gcab = gamma_fixed(&cab, p, cc);
158    let gca = gamma_fixed(&c.sub(a, p, rm()), p, cc);
159    let gcb = gamma_fixed(&c.sub(b, p, rm()), p, cc);
160    if gc.is_nan() || gcab.is_nan() || gca.is_nan() || gcb.is_nan() {
161        return nan_pair(Error::InvalidArgument);
162    }
163    gc.mul(&gcab, p, rm()).div(&gca.mul(&gcb, p, rm()), p, rm())
164}
165
166/// \({}_2F_1(1,1;2;z)=-\ln(1-z)/z\).
167fn is_112(a: &ExactComplex, b: &ExactComplex, c: &ExactComplex, p: usize) -> bool {
168    let one = ExactComplex::one(p);
169    let two = ExactComplex::from_real(ExactNum::from_u8(2, p), p);
170    c_eq(a, &one) && c_eq(b, &one) && c_eq(c, &two)
171}
172
173fn f112(z: &ExactComplex, p: usize, cc: &mut Consts) -> ExactComplex {
174    let one = ExactComplex::one(p);
175    let omz = one.sub(z, p, rm());
176    neg_c(&omz.ln(p, rm(), cc)).div(z, p, rm())
177}
178
179fn pfaff(
180    a: &ExactComplex,
181    b: &ExactComplex,
182    c: &ExactComplex,
183    z: &ExactComplex,
184    p: usize,
185    cc: &mut Consts,
186    left: u32,
187) -> ExactComplex {
188    let one = ExactComplex::one(p);
189    let zm1 = z.sub(&one, p, rm());
190    if is_c_zero(&zm1) {
191        return nan_pair(Error::InvalidArgument);
192    }
193    let w = z.div(&zm1, p, rm());
194    let pref = one.sub(z, p, rm()).pow(&neg_c(a), p, rm(), cc);
195    let cb = c.sub(b, p, rm());
196    pref.mul(&hypergeom_at(a, &cb, c, &w, p, cc, left), p, rm())
197}
198
199/// A&S 15.3.6: argument \(1-z\).
200fn transform_1mz(
201    a: &ExactComplex,
202    b: &ExactComplex,
203    c: &ExactComplex,
204    z: &ExactComplex,
205    p: usize,
206    cc: &mut Consts,
207    left: u32,
208) -> ExactComplex {
209    let one = ExactComplex::one(p);
210    let omz = one.sub(z, p, rm());
211    let cab = c.sub(a, p, rm()).sub(b, p, rm());
212    if is_nonpos_integer(&cab) || is_nonpos_integer(&neg_c(&cab)) {
213        return nan_pair(Error::InvalidArgument);
214    }
215    let gc = gamma_fixed(c, p, cc);
216    let gcab = gamma_fixed(&cab, p, cc);
217    let gca = gamma_fixed(&c.sub(a, p, rm()), p, cc);
218    let gcb = gamma_fixed(&c.sub(b, p, rm()), p, cc);
219    let abc = a.add(b, p, rm()).sub(c, p, rm());
220    let gabc = gamma_fixed(&abc, p, cc);
221    let ga = gamma_fixed(a, p, cc);
222    let gb = gamma_fixed(b, p, cc);
223    if gc.is_nan()
224        || gcab.is_nan()
225        || gca.is_nan()
226        || gcb.is_nan()
227        || gabc.is_nan()
228        || ga.is_nan()
229        || gb.is_nan()
230    {
231        return nan_pair(Error::InvalidArgument);
232    }
233    let a_pref = gc.mul(&gcab, p, rm()).div(&gca.mul(&gcb, p, rm()), p, rm());
234    let b_pref = gc.mul(&gabc, p, rm()).div(&ga.mul(&gb, p, rm()), p, rm());
235    let cap1 = a.add(b, p, rm()).sub(c, p, rm()).add(&one, p, rm());
236    let t1 = a_pref.mul(&hypergeom_at(a, b, &cap1, &omz, p, cc, left), p, rm());
237    let f2 = hypergeom_at(
238        &c.sub(a, p, rm()),
239        &c.sub(b, p, rm()),
240        &cab.add(&one, p, rm()),
241        &omz,
242        p,
243        cc,
244        left,
245    );
246    let t2 = b_pref
247        .mul(&omz.pow(&cab, p, rm(), cc), p, rm())
248        .mul(&f2, p, rm());
249    t1.add(&t2, p, rm())
250}
251
252/// A&S 15.3.7: argument \(1/z\). Degenerate when \(a=b\).
253fn transform_inv(
254    a: &ExactComplex,
255    b: &ExactComplex,
256    c: &ExactComplex,
257    z: &ExactComplex,
258    p: usize,
259    cc: &mut Consts,
260    left: u32,
261) -> ExactComplex {
262    if c_eq(a, b) {
263        return nan_pair(Error::InvalidArgument);
264    }
265    let one = ExactComplex::one(p);
266    let inv = one.div(z, p, rm());
267    let mz = neg_c(z);
268    let t1 = inv_term(a, b, c, &mz, &inv, p, cc, left);
269    let t2 = inv_term(b, a, c, &mz, &inv, p, cc, left);
270    t1.add(&t2, p, rm())
271}
272
273fn inv_term(
274    a: &ExactComplex,
275    b: &ExactComplex,
276    c: &ExactComplex,
277    mz: &ExactComplex,
278    inv: &ExactComplex,
279    p: usize,
280    cc: &mut Consts,
281    left: u32,
282) -> ExactComplex {
283    let one = ExactComplex::one(p);
284    let gc = gamma_fixed(c, p, cc);
285    let gbma = gamma_fixed(&b.sub(a, p, rm()), p, cc);
286    let gb = gamma_fixed(b, p, cc);
287    let gcma = gamma_fixed(&c.sub(a, p, rm()), p, cc);
288    if gc.is_nan() || gbma.is_nan() || gb.is_nan() || gcma.is_nan() {
289        return nan_pair(Error::InvalidArgument);
290    }
291    let pref = gc.mul(&gbma, p, rm()).div(&gb.mul(&gcma, p, rm()), p, rm());
292    let pow = mz.pow(&neg_c(a), p, rm(), cc);
293    let ac1 = a.sub(c, p, rm()).add(&one, p, rm());
294    let ab1 = a.sub(b, p, rm()).add(&one, p, rm());
295    pref.mul(&pow, p, rm())
296        .mul(&hypergeom_at(a, &ac1, &ab1, inv, p, cc, left), p, rm())
297}
298
299fn hypergeom_at(
300    a: &ExactComplex,
301    b: &ExactComplex,
302    c: &ExactComplex,
303    z: &ExactComplex,
304    p: usize,
305    cc: &mut Consts,
306    left: u32,
307) -> ExactComplex {
308    if any_nan(&[a, b, c, z]) {
309        return nan_pair(Error::InvalidArgument);
310    }
311    if is_c_zero(z) || is_c_zero(a) || is_c_zero(b) {
312        return ExactComplex::one(p);
313    }
314    if pole_c(c, a, b, p) {
315        return nan_pair(Error::InvalidArgument);
316    }
317    if is_c_one(z, p) {
318        return kummer_z_one(a, b, c, p, cc);
319    }
320    if is_112(a, b, c, p) {
321        return f112(z, p, cc);
322    }
323    if terminating_neg_int(a, p) || terminating_neg_int(b, p) {
324        return hypergeom_series(a, b, c, z, p);
325    }
326    if abs_lt_one(z, p) {
327        return hypergeom_series(a, b, c, z, p);
328    }
329    if left == 0 {
330        return nan_pair(Error::InvalidArgument);
331    }
332    let next = left - 1;
333    let one = ExactComplex::one(p);
334    let omz = one.sub(z, p, rm());
335    let zm1 = z.sub(&one, p, rm());
336    if re_lt_half(z, p) && !is_c_zero(&zm1) {
337        return pfaff(a, b, c, z, p, cc, next);
338    }
339    if abs_lt_one(&omz, p) && !is_nonpos_integer(&c.sub(a, p, rm()).sub(b, p, rm())) {
340        return transform_1mz(a, b, c, z, p, cc, next);
341    }
342    let inv = one.div(z, p, rm());
343    if abs_lt_one(&inv, p) && !c_eq(a, b) {
344        return transform_inv(a, b, c, z, p, cc, next);
345    }
346    if !is_c_zero(&zm1) {
347        let w = z.div(&zm1, p, rm());
348        if abs_lt_one(&w, p) {
349            return pfaff(a, b, c, z, p, cc, next);
350        }
351    }
352    nan_pair(Error::InvalidArgument)
353}
354
355impl ExactComplex {
356    /// Gaussian \({}_2F_1(a=\mathrm{self},b;c;z)\) in \(\mathbb{C}\).
357    ///
358    /// Series when \(\lvert z\rvert<1\); Pfaff when \(\mathrm{Re}(z)<1/2\);
359    /// Euler / \(1-z\) and \(1/z\) linear transforms otherwise. Kummer at \(z=1\)
360    /// when \(\mathrm{Re}(c-a-b)>0\). Cut on \([1,+\infty)\) in \(z\) (principal
361    /// value from above). Non-positive integer \(c\) (uncanceled) → NaN.
362    ///
363    /// # Precision
364    ///
365    /// - Algorithm: series / Euler / Pfaff / Kummer. Caps `HYPERGEOM_SERIES_MAX_TERMS = 10_000`, `HYPERGEOM_TRANSFORM_MAX = 8`.
366    /// - Bound: Ziv on each part (`MAX_PREC_RETRY`).
367    /// - MPFR oracle: no.
368    pub fn hypergeom_2f1(
369        &self,
370        b: &Self,
371        c: &Self,
372        z: &Self,
373        p: usize,
374        rm: RoundingMode,
375        cc: &mut Consts,
376    ) -> Self {
377        if any_nan(&[self, b, c, z]) {
378            return nan_pair(Error::InvalidArgument);
379        }
380        let dest = round_p(p);
381        if is_c_zero(z) || is_c_zero(self) || is_c_zero(b) {
382            let mut one = ExactComplex::one(dest);
383            one.set_inexact(self.inexact() | b.inexact() | c.inexact() | z.inexact());
384            return one;
385        }
386        if pole_c(c, self, b, dest) {
387            return nan_pair(Error::InvalidArgument);
388        }
389        ziv_complex(dest, rm, |pw| {
390            hypergeom_at(self, b, c, z, pw, cc, HYPERGEOM_TRANSFORM_MAX)
391        })
392    }
393}
394
395#[cfg(test)]
396mod tests {
397    use super::*;
398
399    fn near(a: &ExactNum, b: &ExactNum, p: usize, slack: i32) -> bool {
400        let d = a.sub(b, p, RoundingMode::None).abs();
401        d.is_zero() || d.exponent().unwrap_or(0) < -((p as i32) - slack)
402    }
403
404    fn tiny(x: &ExactNum, p: usize) -> bool {
405        x.is_zero() || x.exponent().unwrap_or(0) < -((p as i32) / 4)
406    }
407
408    #[test]
409    fn complex_hypergeom_2f1_golds() {
410        let p = 256;
411        let r = RoundingMode::ToEven;
412        let mut cc = Consts::new().unwrap();
413        let one = ExactComplex::one(p);
414        let two = ExactComplex::from_real(ExactNum::from_u8(2, p), p);
415        let half = ExactNum::from_u8(1, p).div(&ExactNum::from_u8(2, p), p, r);
416        let zh = ExactComplex::from_real(half.clone(), p);
417        let a_h = zh.clone();
418        let z0 = ExactComplex::zero(p);
419
420        let f1 = one.hypergeom_2f1(&one, &two, &zh, p, r, &mut cc);
421        let ln2 = cc.ln_2(p, r);
422        let want = ln2.mul(&ExactNum::from_u8(2, p), p, r);
423        assert!(near(f1.re(), &want, p, 40), "2F1(1,1;2;1/2)=2ln2");
424        assert!(tiny(f1.im(), p));
425
426        let fk = a_h.hypergeom_2f1(&a_h, &one, &zh, p, r, &mut cc);
427        let k = half.elliptic_k(p, r, &mut cc);
428        let pi = cc.pi(p, r);
429        let want_k = k.mul(&ExactNum::from_u8(2, p), p, r).div(&pi, p, r);
430        assert!(near(fk.re(), &want_k, p, 40), "2F1(1/2,1/2;1;1/2)=2K/π");
431        assert!(tiny(fk.im(), p));
432
433        let a = ExactComplex::new(
434            ExactNum::from_u8(2, p).div(&ExactNum::from_u8(3, p), p, r),
435            ExactNum::from_u8(1, p).div(&ExactNum::from_u8(4, p), p, r),
436        );
437        let b = ExactComplex::new(
438            ExactNum::from_u8(1, p).div(&ExactNum::from_u8(5, p), p, r),
439            ExactNum::from_u8(1, p).div(&ExactNum::from_u8(7, p), p, r),
440        );
441        let c = ExactComplex::new(
442            ExactNum::from_u8(3, p).div(&ExactNum::from_u8(2, p), p, r),
443            ExactNum::from_u8(1, p).div(&ExactNum::from_u8(9, p), p, r),
444        );
445        let f0 = a.hypergeom_2f1(&b, &c, &z0, p, r, &mut cc);
446        assert!(
447            near(f0.re(), &ExactNum::from_u8(1, p), p, 8),
448            "2F1(*,*,*;0)=1"
449        );
450        assert!(tiny(f0.im(), p));
451
452        let z = ExactComplex::new(
453            ExactNum::from_u8(3, p).div(&ExactNum::from_u8(10, p), p, r),
454            ExactNum::from_u8(1, p).div(&ExactNum::from_u8(10, p), p, r),
455        );
456        let lhs = a.hypergeom_2f1(&b, &c, &z, p, r, &mut cc);
457        let cab = c.sub(&a, p, r).sub(&b, p, r);
458        let pref = ExactComplex::one(p).sub(&z, p, r).pow(&cab, p, r, &mut cc);
459        let rhs = pref.mul(
460            &c.sub(&a, p, r)
461                .hypergeom_2f1(&c.sub(&b, p, r), &c, &z, p, r, &mut cc),
462            p,
463            r,
464        );
465        assert!(near(lhs.re(), rhs.re(), p, 30), "Euler re");
466        assert!(near(lhs.im(), rhs.im(), p, 30), "Euler im");
467
468        let c0 = ExactComplex::zero(p);
469        let bad = one.hypergeom_2f1(&one, &c0, &zh, p, r, &mut cc);
470        assert!(bad.is_nan(), "c=0 → NaN");
471
472        let two_r = ExactNum::from_u8(2, p);
473        let eps = ExactNum::from_u8(1, p).ldexp(-40, p, RoundingMode::None);
474        let above = ExactComplex::new(two_r.clone(), eps.clone());
475        let below = ExactComplex::new(two_r, eps.neg());
476        let fa = one.hypergeom_2f1(&one, &two, &above, p, r, &mut cc);
477        let fb = one.hypergeom_2f1(&one, &two, &below, p, r, &mut cc);
478        assert!(!fa.is_nan() && !fb.is_nan(), "cut defined");
479        assert!(near(fa.re(), fb.re(), p, 20), "cut Re");
480        assert_ne!(fa.im().cmp(fb.im()), Some(0), "cut Im differs");
481        assert!(near(&fa.im().abs(), &fb.im().abs(), p, 20), "cut |Im|");
482    }
483}
zenith_float_num/complex_hypergeom.rs

zenith_float_num/
complex_hypergeom.rs
