ellip 1.1.0

Elliptic integrals for Rust
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

/*
 * Ellip is licensed under The 3-Clause BSD, see LICENSE.
 * Copyright 2025 Sira Pornsiriprasert <code@psira.me>
 * This code is modified from Boost Math, see LICENSE in this directory.
 */

// Original header from Boost Math
//  Copyright (c) 2006 Xiaogang Zhang
//  Copyright (c) 2006 John Maddock
//  Copyright (c) 2024 Matt Borland
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0.

use num_traits::Float;

use crate::{carlson::elliprd_unchecked, crate_util::check, StrErr};

/// Computes [complete elliptic integral of Legendre's type](https://dlmf.nist.gov/19.2.E8).
/// ```text
///           π/2
///          ⌠       sin²θ dθ
/// D(m)  =  │  _________________
///          │     _____________
///          ⌡   \╱ 1 - m sin²θ
///         0
/// ```
///
/// ## Parameters
/// - m: elliptic parameter. m ∈ ℝ, m ≤ 1.
///
/// The elliptic modulus (k) is also frequently used instead of the parameter (m), where k² = m.
///
/// ## Domain
/// - Returns error if m > 1.
///
/// ## Graph
/// ![Complete Elliptic Integral of Legendre's Type](https://github.com/p-sira/ellip/blob/main/figures/ellipd.svg?raw=true)
///
/// [Interactive Plot](https://p-sira.github.io/ellippy/_static/figures/ellipd.html)
///
/// ## Special Cases
/// - D(0) = π/4
/// - D(1) = ∞
/// - D(-∞) = 0
///
/// # Related Functions
/// - [ellipd](crate::ellipd)(m) = ([ellipk](crate::ellipk)(m) - [ellipe](crate::ellipe)(m)) / m
/// - [ellipd](crate::ellipd)(m) = [elliprd](crate::elliprd)(0, 1 - m, 1) / 3
/// - [ellipdinc](crate::ellipdinc)(π/2, m) = [ellipd](crate::ellipd)(m)
///
/// # Examples
/// ```
/// use ellip::{ellipd, util::assert_close};
///
/// assert_close(ellipd(0.5).unwrap(), 1.0068615925073927, 1e-15);
/// ```
///
/// # References
/// - Maddock, John, Paul Bristow, Hubert Holin, and Xiaogang Zhang. “Boost Math Library: Special Functions - Elliptic Integrals.” Accessed April 17, 2025. <https://www.boost.org/doc/libs/1_88_0/libs/math/doc/html/math_toolkit/ellint.html>.
/// - Carlson, B. C. “DLMF: Chapter 19 Elliptic Integrals.” Accessed February 19, 2025. <https://dlmf.nist.gov/19>.
#[numeric_literals::replace_float_literals(T::from(literal).unwrap())]
pub fn ellipd<T: Float>(m: T) -> Result<T, StrErr> {
    if m >= 1.0 {
        // D evaluates to inf at m=1.
        if m == 1.0 {
            return Ok(inf!());
        }
        return Err("ellipd: m must not be greater than 1.");
    }

    if m.abs() <= epsilon!() {
        return Ok(pi!() / 4.0);
    }

    // m -> -inf, sqrt(1+m sin²θ) -> m sinθ, θ=[0,pi/2], given phi = pi/2
    // then D(m) -> 1/sqrt(-m)
    if m <= min_val!() {
        return Ok(1.0 / (-m).sqrt());
    }

    let ans = elliprd_unchecked(0.0, 1.0 - m, 1.0) / 3.0;
    #[cfg(feature = "test_force_fail")]
    let ans = nan!();

    if ans.is_finite() {
        return Ok(ans);
    }
    check!(@nan, ellipd, [m]);
    Err("ellipd: Unexpected error.")
}

#[cfg(not(feature = "test_force_fail"))]
#[cfg(test)]
mod tests {
    use super::*;
    use crate::{compare_test_data_boost, compare_test_data_wolfram};

    #[test]
    fn test_ellipd() {
        compare_test_data_boost!("ellipd_data.txt", ellipd, 1, 2.9e-16);
    }

    #[test]
    fn test_ellipd_wolfram() {
        compare_test_data_wolfram!("ellipd_data.csv", ellipd, 1, 8e-16);
        compare_test_data_wolfram!("ellipd_neg.csv", ellipd, 1, 8e-16);
    }

    #[test]
    fn test_ellipd_special_cases() {
        use std::f64::{consts::FRAC_PI_4, INFINITY, MAX, NAN, NEG_INFINITY};
        // m = 0: D(0) = pi/4
        assert_eq!(ellipd(0.0).unwrap(), FRAC_PI_4);
        // m = 1: D(1) = inf
        assert_eq!(ellipd(1.0).unwrap(), INFINITY);
        // m < 0: should be valid
        assert!(ellipd(-1.0).unwrap().is_finite());
        // m > 1: should return Err
        assert_eq!(ellipd(1.1), Err("ellipd: m must not be greater than 1."));
        // m = NaN: should return Err
        assert_eq!(ellipd(NAN), Err("ellipd: Arguments cannot be NAN."));
        // m = inf: should return Err
        assert_eq!(
            ellipd(INFINITY),
            Err("ellipd: m must not be greater than 1.")
        );
        // m -> -inf: D(m) = 1/sqrt(-m)
        assert_eq!(ellipd(-MAX).unwrap(), 1.0 / MAX.sqrt());
        // m = -inf: D(-inf) = 0
        assert_eq!(ellipd(NEG_INFINITY).unwrap(), 0.0);
    }
}

#[cfg(feature = "test_force_fail")]
crate::test_force_unreachable! {
    assert_eq!(ellipd(0.5), Err("ellipd: Unexpected error."));
}