zenith_float_num/ops/
jacobi.rs

1//! Jacobi elliptic functions. Parameter \(m=k^2\). Real \(m\in[0,1]\).
2
3use alloc::vec::Vec;
4
5use crate::common::util::bump_prec_retry;
6use crate::common::util::round_p;
7use crate::defs::Error;
8use crate::defs::Exponent;
9use crate::defs::RoundingMode;
10use crate::defs::Sign;
11use crate::num::ExactNumNumber;
12use crate::ops::consts::Consts;
13use crate::WORD_BIT_SIZE;
14
15/// AGM / descending Landen steps for `am`, `sn`, `cn`, `dn`.
16pub const JACOBI_AGM_MAX: u32 = 128;
17
18impl ExactNumNumber {
19    /// Jacobi amplitude \(\operatorname{am}(u|m)\). Parameter \(m=k^2\in[0,1]\).
20    pub fn jacobi_am(
21        &self,
22        m: &Self,
23        p: usize,
24        rm: RoundingMode,
25        cc: &mut Consts,
26    ) -> Result<Self, Error> {
27        self.jacobi_ziv(m, p, rm, cc, |am, _sn, _cn, _dn, _px| am.clone())
28    }
29
30    /// \(\operatorname{sn}(u|m)=\sin\operatorname{am}(u|m)\).
31    pub fn jacobi_sn(
32        &self,
33        m: &Self,
34        p: usize,
35        rm: RoundingMode,
36        cc: &mut Consts,
37    ) -> Result<Self, Error> {
38        self.jacobi_ziv(m, p, rm, cc, |_am, sn, _cn, _dn, _px| sn.clone())
39    }
40
41    /// \(\operatorname{cn}(u|m)=\cos\operatorname{am}(u|m)\).
42    pub fn jacobi_cn(
43        &self,
44        m: &Self,
45        p: usize,
46        rm: RoundingMode,
47        cc: &mut Consts,
48    ) -> Result<Self, Error> {
49        self.jacobi_ziv(m, p, rm, cc, |_am, _sn, cn, _dn, _px| cn.clone())
50    }
51
52    /// \(\operatorname{dn}(u|m)=\sqrt{1-m\,\operatorname{sn}^2(u|m)}\).
53    pub fn jacobi_dn(
54        &self,
55        m: &Self,
56        p: usize,
57        rm: RoundingMode,
58        cc: &mut Consts,
59    ) -> Result<Self, Error> {
60        self.jacobi_ziv(m, p, rm, cc, |_am, _sn, _cn, dn, _px| dn.clone())
61    }
62
63    /// \(\operatorname{cd}(u|m)=\operatorname{cn}(u|m)/\operatorname{dn}(u|m)\).
64    pub fn jacobi_cd(
65        &self,
66        m: &Self,
67        p: usize,
68        rm: RoundingMode,
69        cc: &mut Consts,
70    ) -> Result<Self, Error> {
71        self.jacobi_ziv(m, p, rm, cc, |_am, _sn, cn, dn, px| div_nz(cn, dn, px))
72    }
73
74    /// \(\operatorname{ns}(u|m)=1/\operatorname{sn}(u|m)\).
75    pub fn jacobi_ns(
76        &self,
77        m: &Self,
78        p: usize,
79        rm: RoundingMode,
80        cc: &mut Consts,
81    ) -> Result<Self, Error> {
82        self.jacobi_ziv(m, p, rm, cc, |_am, sn, _cn, _dn, px| rec_nz(sn, px))
83    }
84
85    /// \(\operatorname{nc}(u|m)=1/\operatorname{cn}(u|m)\).
86    pub fn jacobi_nc(
87        &self,
88        m: &Self,
89        p: usize,
90        rm: RoundingMode,
91        cc: &mut Consts,
92    ) -> Result<Self, Error> {
93        self.jacobi_ziv(m, p, rm, cc, |_am, _sn, cn, _dn, px| rec_nz(cn, px))
94    }
95
96    /// \(\operatorname{nd}(u|m)=1/\operatorname{dn}(u|m)\).
97    pub fn jacobi_nd(
98        &self,
99        m: &Self,
100        p: usize,
101        rm: RoundingMode,
102        cc: &mut Consts,
103    ) -> Result<Self, Error> {
104        self.jacobi_ziv(m, p, rm, cc, |_am, _sn, _cn, dn, px| rec_nz(dn, px))
105    }
106
107    /// \(\operatorname{sc}(u|m)=\operatorname{sn}/\operatorname{cn}\).
108    pub fn jacobi_sc(
109        &self,
110        m: &Self,
111        p: usize,
112        rm: RoundingMode,
113        cc: &mut Consts,
114    ) -> Result<Self, Error> {
115        self.jacobi_ziv(m, p, rm, cc, |_am, sn, cn, _dn, px| div_nz(sn, cn, px))
116    }
117
118    /// \(\operatorname{sd}(u|m)=\operatorname{sn}/\operatorname{dn}\).
119    pub fn jacobi_sd(
120        &self,
121        m: &Self,
122        p: usize,
123        rm: RoundingMode,
124        cc: &mut Consts,
125    ) -> Result<Self, Error> {
126        self.jacobi_ziv(m, p, rm, cc, |_am, sn, _cn, dn, px| div_nz(sn, dn, px))
127    }
128
129    /// \(\operatorname{cs}(u|m)=\operatorname{cn}/\operatorname{sn}\).
130    pub fn jacobi_cs(
131        &self,
132        m: &Self,
133        p: usize,
134        rm: RoundingMode,
135        cc: &mut Consts,
136    ) -> Result<Self, Error> {
137        self.jacobi_ziv(m, p, rm, cc, |_am, sn, cn, _dn, px| div_nz(cn, sn, px))
138    }
139
140    /// \(\operatorname{ds}(u|m)=\operatorname{dn}/\operatorname{sn}\).
141    pub fn jacobi_ds(
142        &self,
143        m: &Self,
144        p: usize,
145        rm: RoundingMode,
146        cc: &mut Consts,
147    ) -> Result<Self, Error> {
148        self.jacobi_ziv(m, p, rm, cc, |_am, sn, _cn, dn, px| div_nz(dn, sn, px))
149    }
150
151    /// \(\operatorname{dc}(u|m)=\operatorname{dn}/\operatorname{cn}\).
152    pub fn jacobi_dc(
153        &self,
154        m: &Self,
155        p: usize,
156        rm: RoundingMode,
157        cc: &mut Consts,
158    ) -> Result<Self, Error> {
159        self.jacobi_ziv(m, p, rm, cc, |_am, _sn, cn, dn, px| div_nz(dn, cn, px))
160    }
161
162    fn jacobi_ziv<F>(
163        &self,
164        m: &Self,
165        p: usize,
166        rm: RoundingMode,
167        cc: &mut Consts,
168        pick: F,
169    ) -> Result<Self, Error>
170    where
171        F: Fn(&Self, &Self, &Self, &Self, usize) -> Result<Self, Error>,
172    {
173        let p = round_p(p);
174        Self::p_assertion(p)?;
175        jacobi_m_in_unit(m, p)?;
176        let mut p_inc = WORD_BIT_SIZE;
177        let mut p_wrk = p
178            .max(self.mantissa_max_bit_len())
179            .max(m.mantissa_max_bit_len())
180            + p_inc;
181        loop {
182            let p_x = p_wrk + WORD_BIT_SIZE * 2;
183            let (am, sn, cn, dn) = self.jacobi_sncndn_at(m, p_x, cc)?;
184            let mut ret = pick(&am, &sn, &cn, &dn, p_x)?;
185            if ret.try_set_precision(p, rm, p_wrk)? {
186                ret.set_inexact(ret.inexact() | self.inexact() | m.inexact());
187                return Ok(ret);
188            }
189            bump_prec_retry(&mut p_wrk, &mut p_inc, p)?;
190        }
191    }
192
193    fn jacobi_sncndn_at(
194        &self,
195        m: &Self,
196        p: usize,
197        cc: &mut Consts,
198    ) -> Result<(Self, Self, Self, Self), Error> {
199        let zero = Self::new2(p, Sign::Pos, false)?;
200        let one = Self::from_word(1, p)?;
201        if self.is_zero() {
202            return Ok((zero.clone()?, zero, one.clone()?, one));
203        }
204        if m.is_zero() {
205            let am = self.clone()?;
206            let sn = self.sin(p, RoundingMode::None, cc)?;
207            let cn = self.cos(p, RoundingMode::None, cc)?;
208            return Ok((am, sn, cn, one));
209        }
210        if m.cmp(&one) == 0 {
211            let sn = self.tanh(p, RoundingMode::None, cc)?;
212            let ch = self.cosh(p, RoundingMode::None, cc)?;
213            let cn = one.div(&ch, p, RoundingMode::None)?;
214            let sh = self.sinh(p, RoundingMode::None, cc)?;
215            let am = sh.atan(p, RoundingMode::None, cc)?;
216            return Ok((am, sn, cn.clone()?, cn));
217        }
218        let u = self.jacobi_reduce_period(m, p, cc)?;
219        jacobi_agm(&u, m, p, cc)
220    }
221
222    fn jacobi_reduce_period(&self, m: &Self, p: usize, cc: &mut Consts) -> Result<Self, Error> {
223        // \(K(m)\ge\pi/2\), so \(2K\ge\pi\). Skip \(K\) (and a nested Ziv) when \(\lvert u\rvert<\pi\).
224        let pi = cc.pi_num(p, RoundingMode::None)?;
225        if self.abs()?.cmp(&pi) < 0 {
226            return self.clone();
227        }
228        let k = m.elliptic_k_at(p, cc)?;
229        let two = Self::from_word(2, p)?;
230        let four = Self::from_word(4, p)?;
231        let two_k = two.mul(&k, p, RoundingMode::None)?;
232        if self.abs()?.cmp(&two_k) <= 0 {
233            return self.clone();
234        }
235        let four_k = four.mul(&k, p, RoundingMode::None)?;
236        let mut r = self.rem(&four_k)?;
237        if r.cmp(&two_k) > 0 {
238            r = r.sub(&four_k, p, RoundingMode::None)?;
239        } else if r.is_negative() && r.abs()?.cmp(&two_k) > 0 {
240            r = r.add(&four_k, p, RoundingMode::None)?;
241        }
242        Ok(r)
243    }
244}
245
246fn div_nz(n: &ExactNumNumber, d: &ExactNumNumber, p: usize) -> Result<ExactNumNumber, Error> {
247    if d.is_zero() {
248        return Err(Error::InvalidArgument);
249    }
250    n.div(d, p, RoundingMode::None)
251}
252
253fn rec_nz(d: &ExactNumNumber, p: usize) -> Result<ExactNumNumber, Error> {
254    if d.is_zero() {
255        return Err(Error::InvalidArgument);
256    }
257    ExactNumNumber::from_word(1, p)?.div(d, p, RoundingMode::None)
258}
259
260fn jacobi_m_in_unit(m: &ExactNumNumber, p: usize) -> Result<(), Error> {
261    if m.is_negative() {
262        return Err(Error::InvalidArgument);
263    }
264    let one = ExactNumNumber::from_word(1, p)?;
265    if m.cmp(&one) > 0 {
266        return Err(Error::InvalidArgument);
267    }
268    Ok(())
269}
270
271fn agm_small(c: &ExactNumNumber, a: &ExactNumNumber, p: usize) -> Result<bool, Error> {
272    if c.is_zero() {
273        return Ok(true);
274    }
275    if a.is_zero() {
276        return Ok(false);
277    }
278    let rel = c.abs()?.div(&a.abs()?, p, RoundingMode::None)?;
279    let thresh = -((p as i32) / 2 + 16);
280    Ok(rel.exponent() < thresh)
281}
282
283fn clamp_unit(x: &ExactNumNumber, p: usize) -> Result<ExactNumNumber, Error> {
284    let one = ExactNumNumber::from_word(1, p)?;
285    let ax = x.abs()?;
286    if ax.cmp(&one) <= 0 {
287        return x.clone();
288    }
289    let mut y = one;
290    y.set_sign(x.sign());
291    Ok(y)
292}
293
294fn jacobi_agm(
295    u: &ExactNumNumber,
296    m: &ExactNumNumber,
297    p: usize,
298    cc: &mut Consts,
299) -> Result<
300    (
301        ExactNumNumber,
302        ExactNumNumber,
303        ExactNumNumber,
304        ExactNumNumber,
305    ),
306    Error,
307> {
308    let one = ExactNumNumber::from_word(1, p)?;
309    let two = ExactNumNumber::from_word(2, p)?;
310    let om = one.sub(m, p, RoundingMode::None)?;
311    if om.is_negative() {
312        return Err(Error::InvalidArgument);
313    }
314    let mut b = om.sqrt(p, RoundingMode::None)?;
315    let mut aa = Vec::new();
316    let mut cc_seq = Vec::new();
317    aa.push(one.clone()?);
318    cc_seq.push(m.sqrt(p, RoundingMode::None)?);
319    let mut n = 0u32;
320    while n < JACOBI_AGM_MAX {
321        if agm_small(&cc_seq[n as usize], &aa[n as usize], p)? {
322            break;
323        }
324        let ai = aa[n as usize].clone()?;
325        let c_next = ai
326            .sub(&b, p, RoundingMode::None)?
327            .div(&two, p, RoundingMode::None)?;
328        let t = ai
329            .mul(&b, p, RoundingMode::None)?
330            .sqrt(p, RoundingMode::None)?;
331        let a_next = ai
332            .add(&b, p, RoundingMode::None)?
333            .div(&two, p, RoundingMode::None)?;
334        b = t;
335        aa.push(a_next);
336        cc_seq.push(c_next);
337        n += 1;
338    }
339    if n == JACOBI_AGM_MAX && !agm_small(&cc_seq[n as usize], &aa[n as usize], p)? {
340        return Err(Error::InvalidArgument);
341    }
342    let last = aa.len() - 1;
343    let shift = Exponent::try_from(last).map_err(|_| Error::InvalidArgument)?;
344    let twon = one.ldexp(shift, p, RoundingMode::None)?;
345    let mut phi = twon
346        .mul(&aa[last], p, RoundingMode::None)?
347        .mul(u, p, RoundingMode::None)?;
348    let mut idx = last;
349    while idx > 0 {
350        let sphi = phi.sin(p, RoundingMode::None, cc)?;
351        let t = clamp_unit(
352            &cc_seq[idx]
353                .mul(&sphi, p, RoundingMode::None)?
354                .div(&aa[idx], p, RoundingMode::None)?,
355            p,
356        )?;
357        let asint = t.asin(p, RoundingMode::None, cc)?;
358        phi = asint
359            .add(&phi, p, RoundingMode::None)?
360            .div(&two, p, RoundingMode::None)?;
361        idx -= 1;
362    }
363    let sn = phi.sin(p, RoundingMode::None, cc)?;
364    let cn = phi.cos(p, RoundingMode::None, cc)?;
365    let msn2 = m.mul(&sn.mul(&sn, p, RoundingMode::None)?, p, RoundingMode::None)?;
366    let dnarg = one.sub(&msn2, p, RoundingMode::None)?;
367    if dnarg.is_negative() {
368        return Err(Error::InvalidArgument);
369    }
370    let dn = dnarg.sqrt(p, RoundingMode::None)?;
371    Ok((phi, sn, cn, dn))
372}
373
374#[cfg(test)]
375mod tests {
376    use super::*;
377    use crate::defs::RoundingMode;
378
379    fn bits_agree(a: &ExactNumNumber, b: &ExactNumNumber, p: usize, min_bits: i32, label: &str) {
380        let d = a.sub(b, p, RoundingMode::None).unwrap().abs().unwrap();
381        if d.is_zero() {
382            return;
383        }
384        let rel = d.exponent() - a.exponent();
385        assert!(
386            rel < -min_bits,
387            "{label}: relative exponent {rel} (a.exp={}, b.exp={})",
388            a.exponent(),
389            b.exponent()
390        );
391    }
392
393    #[test]
394    fn test_jacobi_identities() {
395        let p = 256;
396        let mut cc = Consts::new().unwrap();
397        let rm = RoundingMode::ToEven;
398        let one = ExactNumNumber::from_word(1, p).unwrap();
399        let zero = ExactNumNumber::new(p).unwrap();
400        let two = ExactNumNumber::from_word(2, p).unwrap();
401        let mut half = one.clone().unwrap();
402        half.div_by_2(RoundingMode::None);
403
404        let sn0 = zero.jacobi_sn(&half, p, rm, &mut cc).unwrap();
405        assert!(sn0.is_zero(), "sn(0)");
406        let cn0 = zero.jacobi_cn(&half, p, rm, &mut cc).unwrap();
407        bits_agree(&cn0, &one, p, 80, "cn(0)");
408        let dn0 = zero.jacobi_dn(&half, p, rm, &mut cc).unwrap();
409        bits_agree(&dn0, &one, p, 80, "dn(0)");
410        let am0 = zero.jacobi_am(&half, p, rm, &mut cc).unwrap();
411        assert!(am0.is_zero(), "am(0)");
412
413        let u = one.clone().unwrap();
414        let sn_m0 = u.jacobi_sn(&zero, p, rm, &mut cc).unwrap();
415        let sinu = u.sin(p, rm, &mut cc).unwrap();
416        bits_agree(&sn_m0, &sinu, p, 80, "sn(u|0)=sin u");
417        let cn_m0 = u.jacobi_cn(&zero, p, rm, &mut cc).unwrap();
418        let cosu = u.cos(p, rm, &mut cc).unwrap();
419        bits_agree(&cn_m0, &cosu, p, 80, "cn(u|0)=cos u");
420        let dn_m0 = u.jacobi_dn(&zero, p, rm, &mut cc).unwrap();
421        bits_agree(&dn_m0, &one, p, 80, "dn(u|0)=1");
422
423        let sn_m1 = u.jacobi_sn(&one, p, rm, &mut cc).unwrap();
424        let tanhu = u.tanh(p, rm, &mut cc).unwrap();
425        bits_agree(&sn_m1, &tanhu, p, 80, "sn(u|1)=tanh u");
426        let cn_m1 = u.jacobi_cn(&one, p, rm, &mut cc).unwrap();
427        let sech = one
428            .div(&u.cosh(p, rm, &mut cc).unwrap(), p, RoundingMode::None)
429            .unwrap();
430        bits_agree(&cn_m1, &sech, p, 80, "cn(u|1)=sech u");
431        let dn_m1 = u.jacobi_dn(&one, p, rm, &mut cc).unwrap();
432        bits_agree(&dn_m1, &sech, p, 80, "dn(u|1)=sech u");
433
434        let sn = u.jacobi_sn(&half, p, rm, &mut cc).unwrap();
435        let cn = u.jacobi_cn(&half, p, rm, &mut cc).unwrap();
436        let dn = u.jacobi_dn(&half, p, rm, &mut cc).unwrap();
437        let s2c2 = sn
438            .mul(&sn, p, RoundingMode::None)
439            .unwrap()
440            .add(
441                &cn.mul(&cn, p, RoundingMode::None).unwrap(),
442                p,
443                RoundingMode::None,
444            )
445            .unwrap();
446        bits_agree(&s2c2, &one, p, 80, "sn²+cn²=1");
447        let d2ms = dn
448            .mul(&dn, p, RoundingMode::None)
449            .unwrap()
450            .add(
451                &half
452                    .mul(
453                        &sn.mul(&sn, p, RoundingMode::None).unwrap(),
454                        p,
455                        RoundingMode::None,
456                    )
457                    .unwrap(),
458                p,
459                RoundingMode::None,
460            )
461            .unwrap();
462        bits_agree(&d2ms, &one, p, 80, "dn²+m sn²=1");
463
464        let k = half.elliptic_k(p, rm, &mut cc).unwrap();
465        let mut kh = k.clone().unwrap();
466        kh.div_by_2(RoundingMode::None);
467        let snkh = kh.jacobi_sn(&half, p, rm, &mut cc).unwrap();
468        let kp = one
469            .sub(&half, p, RoundingMode::None)
470            .unwrap()
471            .sqrt(p, RoundingMode::None)
472            .unwrap();
473        let want = one
474            .div(
475                &one.add(&kp, p, RoundingMode::None)
476                    .unwrap()
477                    .sqrt(p, RoundingMode::None)
478                    .unwrap(),
479                p,
480                RoundingMode::None,
481            )
482            .unwrap();
483        bits_agree(&snkh, &want, p, 40, "sn(K/2|1/2)");
484
485        let four_k = four_word(p).mul(&k, p, RoundingMode::None).unwrap();
486        let u4 = u.add(&four_k, p, RoundingMode::None).unwrap();
487        let sn_per = u4.jacobi_sn(&half, p, rm, &mut cc).unwrap();
488        bits_agree(&sn_per, &sn, p, 40, "sn(u+4K)=sn(u)");
489
490        let mut h = one.clone().unwrap();
491        for _ in 0..20 {
492            h.div_by_2(RoundingMode::None);
493        }
494        let up = u.add(&h, p, RoundingMode::None).unwrap();
495        let um = u.sub(&h, p, RoundingMode::None).unwrap();
496        let snp = up.jacobi_sn(&half, p, rm, &mut cc).unwrap();
497        let snm = um.jacobi_sn(&half, p, rm, &mut cc).unwrap();
498        let fd = snp
499            .sub(&snm, p, RoundingMode::None)
500            .unwrap()
501            .div(
502                &two.mul(&h, p, RoundingMode::None).unwrap(),
503                p,
504                RoundingMode::None,
505            )
506            .unwrap();
507        let deriv = cn.mul(&dn, p, RoundingMode::None).unwrap();
508        bits_agree(&fd, &deriv, p, 20, "d(sn)/du = cn dn");
509
510        let sn_none = u
511            .jacobi_sn(&half, 512, RoundingMode::None, &mut cc)
512            .unwrap();
513        bits_agree(&sn_none, &sn, p, 20, "sn None 512 vs ToEven 256");
514
515        assert!(u.jacobi_sn(&one.neg().unwrap(), p, rm, &mut cc).is_err());
516        assert!(u.jacobi_sn(&two, p, rm, &mut cc).is_err());
517        let cd = u.jacobi_cd(&half, p, rm, &mut cc).unwrap();
518        let want_cd = cn.div(&dn, p, RoundingMode::None).unwrap();
519        bits_agree(&cd, &want_cd, p, 80, "cd=cn/dn");
520        let ns = u.jacobi_ns(&half, p, rm, &mut cc).unwrap();
521        let nssn = ns.mul(&sn, p, RoundingMode::None).unwrap();
522        bits_agree(&nssn, &one, p, 80, "ns·sn=1");
523        assert!(zero.jacobi_ns(&half, p, rm, &mut cc).is_err());
524
525        // Inverse of sn is incomplete F: F(sn(u|m)|m)=u for |u|<K(m).
526        let back = sn.elliptic_f(&half, p, rm, &mut cc).unwrap();
527        bits_agree(&back, &u, p, 40, "F(sn(u|m)|m)=u");
528    }
529
530    fn four_word(p: usize) -> ExactNumNumber {
531        ExactNumNumber::from_word(4, p).unwrap()
532    }
533}
zenith_float_num/ops/jacobi.rs

zenith_float_num/ops/
jacobi.rs
