/// Get a function that returns the derivative of the provided univariate, scalar-valued function.
///
/// The derivative is computed using forward-mode automatic differentiation.
///
/// # Arguments
///
/// * `f` - Univariate, scalar-valued function, $f:\mathbb{R}\to\mathbb{R}$.
/// * `func_name` - Name of the function that will return the derivative of $f(x)$ at any point
///   $x\in\mathbb{R}$.
/// * `param_type` (optional) - Type of the extra runtime parameter `p` that is passed to `f`.
///   Defaults to `[f64]` (implying that `f` accepts `p: &[f64]`).
///
/// # Warning
///
/// `f` cannot be defined as closure. It must be defined as a function.
///
/// # Note
///
/// The function produced by this macro will perform 1 evaluation of $f(x)$ to evaluate its
/// derivative.
///
/// # Examples
///
/// ## Basic Example
///
/// Compute the derivative of
///
/// $$f(x)=x^{3}$$
///
/// at $x=2$, and compare the result to the true result of $f'(2)=12$.
///
/// ```
/// use linalg_traits::Scalar;
/// use numtest::*;
///
/// use numdiff::{get_sderivative, Dual};
///
/// // Define the function, f(x).
/// fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
///     x.powi(3)
/// }
///
/// // Parameter vector (empty for this example).
/// let p = [];
///
/// // Autogenerate the function "df" that can be used to compute the derivative of f(x) at any
/// // point x.
/// get_sderivative!(f, df);
///
/// // Compute the derivative of f(x) at the evaluation point, x = 2.
/// let df_at_2: f64 = df(2.0, &p);
///
/// // Check the accuracy of the derivative.
/// assert_equal_to_decimal!(df_at_2, 12.0, 16);
/// ```
///
/// ## Example Passing Runtime Parameters
///
/// Compute the derivative of a parameterized function
///
/// $$f(x)=ax^{2}+bx+c$$
///
/// where $a$, $b$, and $c$ are runtime parameters. Compare the result against the true derivative
/// of
///
/// $$f'(x)=2ax+b$$
///
/// ```
/// use linalg_traits::Scalar;
/// use numtest::*;
///
/// use numdiff::{get_sderivative, Dual};
///
/// // Define the function, f(x).
/// fn f<S: Scalar>(x: S, p: &[f64]) -> S {
///     let a = S::new(p[0]);
///     let b = S::new(p[1]);
///     let c = S::new(p[2]);
///     a * x.powi(2) + b * x + c
/// }
///
/// // Parameter vector.
/// let a = 2.5;
/// let b = -1.3;
/// let c = 4.7;
/// let p = [a, b, c];
///
/// // Autogenerate the derivative function.
/// get_sderivative!(f, df);
///
/// // True derivative function.
/// let df_true = |x: f64| 2.0 * a * x + b;
///
/// // Compute the derivative at x = 1.0 using both the automatically generated derivative function
/// // and the true derivative function, and compare the results.
/// let df_at_1: f64 = df(1.0, &p);
/// let df_at_1_true: f64 = df_true(1.0);
/// assert_eq!(df_at_1, df_at_1_true);
/// ```
///
/// ## Example Passing Custom Parameter Types
///
/// Use a custom parameter struct instead of `f64` values.
///
/// ```
/// use linalg_traits::Scalar;
/// use numtest::*;
///
/// use numdiff::{get_sderivative, Dual};
///
/// struct Data {
///     a: f64,
///     b: f64,
///     c: f64,
/// }
///
/// // Define the function, f(x).
/// fn f<S: Scalar>(x: S, p: &Data) -> S {
///     let a = S::new(p.a);
///     let b = S::new(p.b);
///     let c = S::new(p.c);
///     a * x.powi(2) + b * x + c
/// }
///
/// // Runtime parameter struct.
/// let p = Data {
///     a: 2.5,
///     b: -1.3,
///     c: 4.7,
/// };
///
/// // Autogenerate the derivative function, telling the macro to expect a runtime parameter of type
/// // &Data.
/// get_sderivative!(f, df, Data);
///
/// // True derivative function.
/// let df_true = |x: f64| 2.0 * p.a * x + p.b;
///
/// // Compute the derivative at x = 1.0 using both the automatically generated derivative function
/// // and the true derivative function, and compare the results.
/// let x0 = 1.0;
/// let df_eval: f64 = df(x0, &p);
/// let df_eval_true: f64 = df_true(x0);
/// assert_equal_to_decimal!(df_eval, df_eval_true, 15);
/// ```
#[macro_export]
macro_rules! get_sderivative {
    ($f:ident, $func_name:ident) => {
        get_sderivative!($f, $func_name, [f64]);
    };
    ($f:ident, $func_name:ident, $param_type:ty) => {
        /// Derivative of a univariate, scalar-valued function `f: ℝ → ℝ`.
        ///
        /// This function is generated for a specific function `f` using the
        /// `numdiff::get_sderivative!` macro.
        ///
        /// # Arguments
        ///
        /// * `x0` - Evaluation point, `x₀ ∈ ℝ`.
        /// * `p` - Extra runtime parameter. This is a parameter (can be of any arbitrary type)
        ///   defined at runtime that the function may depend on but is not differentiated with
        ///   respect to.
        ///
        /// # Returns
        ///
        /// Derivative of `f` with respect to `x`, evaluated at `x = x₀`.
        ///
        /// `(df/dx)|ₓ₌ₓ₀ ∈ ℝ`
        fn $func_name<S: Scalar>(x0: S, p: &$param_type) -> f64 {
            // Step forward in the dual direction.
            let x0 = Dual::new(x0.to_f64().unwrap(), 1.0);

            // Evaluate the function at the dual number.
            let f_x0 = $f(x0, p);

            // Derivative of f with respect to x evaluated at x = x₀.
            f_x0.get_dual()
        }
    };
}

#[cfg(test)]
mod tests {
    use crate::{Dual, test_utils};
    use linalg_traits::Scalar;
    use numtest::*;
    use std::f64::consts::PI;

    #[cfg(feature = "trig")]
    use trig::Trig;

    #[test]
    fn test_product_rule() {
        // f(x) and f'(x).
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(3)
        }
        fn df<S: Scalar>(x: S, _p: &[f64]) -> S {
            S::new(3.0) * x.powi(2)
        }

        // g(x) and g'(x)
        fn g<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.sin()
        }
        fn dg<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.cos()
        }

        // h(x) = f(x)g(x) and h'(x) = f'(x)g(x) + f(x)g'(x).
        fn h<S: Scalar>(x: S, p: &[f64]) -> S {
            f(x, p) * g(x, p)
        }
        fn dh<S: Scalar>(x: S, p: &[f64]) -> S {
            df(x, p) * g(x, p) + f(x, p) * dg(x, p)
        }

        // Parameter vector (empty for this test).
        let p = [];

        // h'(x) obtained via forward-mode automatic differentiation.
        get_sderivative!(h, dh_autodiff);

        // Test autodiff h'(x) against true h'(x).
        assert_eq!(dh_autodiff(-1.5, &p), dh(-1.5, &p));
        assert_eq!(dh_autodiff(1.5, &p), dh(1.5, &p));
    }

    #[test]
    fn test_quotient_rule() {
        // f(x) and f'(x).
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(3)
        }
        fn df<S: Scalar>(x: S, _p: &[f64]) -> S {
            S::new(3.0) * x.powi(2)
        }

        // g(x) and g'(x)
        fn g<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.sin()
        }
        fn dg<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.cos()
        }

        // h(x) = f(x) / g(x) and h'(x) = (g(x)f'(x) - f(x)g'(x)) / (g(x))².
        fn h<S: Scalar>(x: S, p: &[f64]) -> S {
            f(x, p) / g(x, p)
        }
        fn dh<S: Scalar>(x: S, p: &[f64]) -> S {
            (g(x, p) * df(x, p) - f(x, p) * dg(x, p)) / g(x, p).powi(2)
        }

        // Parameter vector (empty for this test).
        let p = [];

        // h'(x) obtained via forward-mode automatic differentiation.
        get_sderivative!(h, dh_autodiff);

        // Test autodiff h'(x) against true h'(x).
        assert_eq!(dh_autodiff(-1.5, &p), dh(-1.5, &p));
        assert_eq!(dh_autodiff(1.5, &p), dh(1.5, &p));
    }

    #[test]
    fn test_chain_rule_one_composition() {
        // f(x) and f'(x).
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(3)
        }
        fn df<S: Scalar>(x: S) -> S {
            S::new(3.0) * x.powi(2)
        }

        // g(x) and g'(x)
        fn g<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.sin()
        }
        fn dg<S: Scalar>(x: S) -> S {
            x.cos()
        }

        // h(x) = g(f(x)) and h'(x) = [g'(f(x))][f'(x)].
        fn h<S: Scalar>(x: S, p: &[f64]) -> S {
            g(f(x, p), p)
        }
        fn dh<S: Scalar>(x: S) -> S {
            dg(f(x, &[])) * df(x)
        }

        // h'(x) obtained via forward-mode automatic differentiation.
        get_sderivative!(h, dh_autodiff);

        // Test autodiff h'(x) against true h'(x).
        assert_eq!(dh_autodiff(-1.5, &[]), dh(-1.5));
        assert_eq!(dh_autodiff(1.5, &[]), dh(1.5));
    }

    #[test]
    fn test_chain_rule_two_compositions() {
        // f(x) and f'(x).
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(3)
        }
        fn df<S: Scalar>(x: S) -> S {
            S::new(3.0) * x.powi(2)
        }

        // g(x) and g'(x)
        fn g<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.sin()
        }
        fn dg<S: Scalar>(x: S) -> S {
            x.cos()
        }

        // h(x) and h'(x).
        fn h<S: Scalar>(x: S, _p: &[f64]) -> S {
            S::new(5.0) / x.powi(2)
        }
        fn dh<S: Scalar>(x: S) -> S {
            S::new(-10.0) / x.powi(3)
        }

        // j(x) = h(g(f(x))) and j'(x) = [h'(g(f(x)))][g'(f(x))][f'(x)].
        fn j<S: Scalar>(x: S, p: &[f64]) -> S {
            h(g(f(x, p), p), p)
        }
        fn dj<S: Scalar>(x: S) -> S {
            dh(g(f(x, &[]), &[])) * dg(f(x, &[])) * df(x)
        }

        // j'(x) obtained via forward-mode automatic differentiation.
        get_sderivative!(j, dj_autodiff);

        // Test autodiff j'(x) against true j'(x).
        assert_equal_to_decimal!(dj_autodiff(1.5, &[]), dj(1.5), 12);
    }

    #[test]
    #[allow(clippy::items_after_statements)]
    fn test_sderivative_polynomial() {
        // Test #1.
        fn f1<S: Scalar>(_x: S, _p: &[f64]) -> S {
            S::one()
        }
        get_sderivative!(f1, df1);
        assert_eq!(df1(2.0, &[]), test_utils::polyi_deriv(0, 2.0));

        // Test #2.
        fn f2<S: Scalar>(x: S, _p: &[f64]) -> S {
            x
        }
        get_sderivative!(f2, df2);
        assert_eq!(df2(2.0, &[]), test_utils::polyi_deriv(1, 2.0));

        // Test #3.
        fn f3<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(2)
        }
        get_sderivative!(f3, df3);
        assert_eq!(df3(2.0, &[]), test_utils::polyi_deriv(2, 2.0));

        // Test #4.
        fn f4<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(3)
        }
        get_sderivative!(f4, df4);
        assert_eq!(df4(2.0, &[]), test_utils::polyi_deriv(3, 2.0));

        // Test #5.
        fn f5<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(4)
        }
        get_sderivative!(f5, df5);
        assert_eq!(df5(2.0, &[]), test_utils::polyi_deriv(4, 2.0));

        // Test #6.
        fn f6<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(5)
        }
        get_sderivative!(f6, df6);
        assert_eq!(df6(2.0, &[]), test_utils::polyi_deriv(5, 2.0));

        // Test #7.
        fn f7<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(6)
        }
        get_sderivative!(f7, df7);
        assert_eq!(df7(2.0, &[]), test_utils::polyi_deriv(6, 2.0));

        // Test #8.
        fn f8<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(7)
        }
        get_sderivative!(f8, df8);
        assert_eq!(df8(2.0, &[]), test_utils::polyi_deriv(7, 2.0));

        // Test #9.
        fn f9<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(-1)
        }
        get_sderivative!(f9, df9);
        assert_eq!(df9(2.0, &[]), test_utils::polyi_deriv(-1, 2.0));

        // Test #10.
        fn f10<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(-2)
        }
        get_sderivative!(f10, df10);
        assert_eq!(df10(2.0, &[]), test_utils::polyi_deriv(-2, 2.0));

        // Test #11.
        fn f11<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(-3)
        }
        get_sderivative!(f11, df11);
        assert_eq!(df11(2.0, &[]), test_utils::polyi_deriv(-3, 2.0));

        // Test #12.
        fn f12<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powi(-7)
        }
        get_sderivative!(f12, df12);
        assert_eq!(df12(2.0, &[]), test_utils::polyi_deriv(-7, 2.0));

        // Test #13.
        fn f13<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powf(S::new(1.0 / 3.0))
        }
        get_sderivative!(f13, df13);
        assert_eq!(df13(2.0, &[]), test_utils::polyf_deriv(1.0 / 3.0, 2.0));

        // Test #14.
        fn f14<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powf(S::new(7.0 / 3.0))
        }
        get_sderivative!(f14, df14);
        assert_equal_to_decimal!(df14(2.0, &[]), test_utils::polyf_deriv(7.0 / 3.0, 2.0), 15);

        // Test #15.
        fn f15<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powf(S::new(-1.0 / 3.0))
        }
        get_sderivative!(f15, df15);
        assert_eq!(df15(2.0, &[]), test_utils::polyf_deriv(-1.0 / 3.0, 2.0));

        // Test #16.
        fn f16<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.powf(S::new(-7.0 / 3.0))
        }
        get_sderivative!(f16, df16);
        assert_eq!(df16(2.0, &[]), test_utils::polyf_deriv(-7.0 / 3.0, 2.0));
    }

    #[test]
    fn test_sderivative_square_root() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.sqrt()
        }
        get_sderivative!(f, df);
        assert_eq!(df(0.5, &[]), test_utils::sqrt_deriv(0.5));
        assert_eq!(df(1.5, &[]), test_utils::sqrt_deriv(1.5));
    }

    #[test]
    fn test_sderivative_exponential() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.exp()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-1.0, &[]), test_utils::exp_deriv(-1.0));
        assert_eq!(df(0.0, &[]), test_utils::exp_deriv(0.0));
        assert_eq!(df(1.0, &[]), test_utils::exp_deriv(1.0));
    }

    #[test]
    fn test_sderivative_power() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            S::new(5.0).powf(x)
        }
        get_sderivative!(f, df);
        assert_eq!(df(-1.0, &[]), test_utils::power_deriv(5.0, -1.0));
        assert_eq!(df(0.0, &[]), test_utils::power_deriv(5.0, 0.0));
        assert_equal_to_decimal!(df(1.0, &[]), test_utils::power_deriv(5.0, 1.0), 14);
    }

    #[test]
    fn test_sderivative_natural_logarithm() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.ln()
        }
        get_sderivative!(f, df);
        assert_eq!(df(0.5, &[]), test_utils::ln_deriv(0.5));
        assert_eq!(df(1.0, &[]), test_utils::ln_deriv(1.0));
        assert_eq!(df(1.5, &[]), test_utils::ln_deriv(1.5));
    }

    #[test]
    fn test_sderivative_base_10_logarithm() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.log10()
        }
        get_sderivative!(f, df);
        assert_eq!(df(0.5, &[]), test_utils::log10_deriv(0.5));
        assert_eq!(df(1.0, &[]), test_utils::log10_deriv(1.0));
        assert_eq!(df(1.5, &[]), test_utils::log10_deriv(1.5));
    }

    #[test]
    fn test_sderivative_sine() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.sin()
        }
        get_sderivative!(f, df);
        assert_eq!(df(0.0, &[]), test_utils::sin_deriv(0.0));
        assert_eq!(df(PI / 4.0, &[]), test_utils::sin_deriv(PI / 4.0));
        assert_eq!(df(PI / 2.0, &[]), test_utils::sin_deriv(PI / 2.0));
        assert_eq!(
            df(3.0 * PI / 4.0, &[]),
            test_utils::sin_deriv(3.0 * PI / 4.0)
        );
        assert_eq!(df(PI, &[]), test_utils::sin_deriv(PI));
        assert_eq!(
            df(5.0 * PI / 4.0, &[]),
            test_utils::sin_deriv(5.0 * PI / 4.0)
        );
        assert_eq!(
            df(3.0 * PI / 2.0, &[]),
            test_utils::sin_deriv(3.0 * PI / 2.0)
        );
        assert_eq!(
            df(7.0 * PI / 4.0, &[]),
            test_utils::sin_deriv(7.0 * PI / 4.0)
        );
        assert_eq!(df(2.0 * PI, &[]), test_utils::sin_deriv(2.0 * PI));
    }

    #[test]
    fn test_sderivative_cosine() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.cos()
        }
        get_sderivative!(f, df);
        assert_eq!(df(0.0, &[]), test_utils::cos_deriv(0.0));
        assert_eq!(df(PI / 4.0, &[]), test_utils::cos_deriv(PI / 4.0));
        assert_eq!(df(PI / 2.0, &[]), test_utils::cos_deriv(PI / 2.0));
        assert_eq!(
            df(3.0 * PI / 4.0, &[]),
            test_utils::cos_deriv(3.0 * PI / 4.0)
        );
        assert_eq!(df(PI, &[]), test_utils::cos_deriv(PI));
        assert_eq!(
            df(5.0 * PI / 4.0, &[]),
            test_utils::cos_deriv(5.0 * PI / 4.0)
        );
        assert_eq!(
            df(3.0 * PI / 2.0, &[]),
            test_utils::cos_deriv(3.0 * PI / 2.0)
        );
        assert_eq!(
            df(7.0 * PI / 4.0, &[]),
            test_utils::cos_deriv(7.0 * PI / 4.0)
        );
        assert_eq!(df(2.0 * PI, &[]), test_utils::cos_deriv(2.0 * PI));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_tangent() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.tan()
        }
        get_sderivative!(f, df);
        assert_eq!(df(0.0, &[]), test_utils::tan_deriv(0.0));
        assert_eq!(df(PI / 4.0, &[]), test_utils::tan_deriv(PI / 4.0));
        assert_eq!(
            df(3.0 * PI / 4.0, &[]),
            test_utils::tan_deriv(3.0 * PI / 4.0)
        );
        assert_eq!(df(PI, &[]), test_utils::tan_deriv(PI));
        assert_eq!(
            df(5.0 * PI / 4.0, &[]),
            test_utils::tan_deriv(5.0 * PI / 4.0)
        );
        assert_eq!(
            df(7.0 * PI / 4.0, &[]),
            test_utils::tan_deriv(7.0 * PI / 4.0)
        );
        assert_eq!(df(2.0 * PI, &[]), test_utils::tan_deriv(2.0 * PI));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_cosecant() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.csc()
        }
        get_sderivative!(f, df);
        assert_equal_to_decimal!(df(PI / 4.0, &[]), test_utils::csc_deriv(PI / 4.0), 15);
        assert_equal_to_decimal!(
            df(3.0 * PI / 4.0, &[]),
            test_utils::csc_deriv(3.0 * PI / 4.0),
            15
        );
        assert_equal_to_decimal!(
            df(5.0 * PI / 4.0, &[]),
            test_utils::csc_deriv(5.0 * PI / 4.0),
            15
        );
        assert_equal_to_decimal!(
            df(7.0 * PI / 4.0, &[]),
            test_utils::csc_deriv(7.0 * PI / 4.0),
            15
        );
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_secant() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.sec()
        }
        get_sderivative!(f, df);
        assert_eq!(df(0.0, &[]), test_utils::sec_deriv(0.0));
        assert_eq!(df(PI / 4.0, &[]), test_utils::sec_deriv(PI / 4.0));
        assert_eq!(
            df(3.0 * PI / 4.0, &[]),
            test_utils::sec_deriv(3.0 * PI / 4.0)
        );
        assert_eq!(df(PI, &[]), test_utils::sec_deriv(PI));
        assert_eq!(
            df(5.0 * PI / 4.0, &[]),
            test_utils::sec_deriv(5.0 * PI / 4.0)
        );
        assert_eq!(
            df(7.0 * PI / 4.0, &[]),
            test_utils::sec_deriv(7.0 * PI / 4.0)
        );
        assert_eq!(df(2.0 * PI, &[]), test_utils::sec_deriv(2.0 * PI));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_cotangent() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.cot()
        }
        get_sderivative!(f, df);
        assert_equal_to_decimal!(df(PI / 4.0, &[]), test_utils::cot_deriv(PI / 4.0), 15);
        assert_equal_to_decimal!(df(PI / 2.0, &[]), test_utils::cot_deriv(PI / 2.0), 16);
        assert_eq!(
            df(3.0 * PI / 4.0, &[]),
            test_utils::cot_deriv(3.0 * PI / 4.0)
        );
        assert_eq!(
            df(5.0 * PI / 4.0, &[]),
            test_utils::cot_deriv(5.0 * PI / 4.0)
        );
        assert_equal_to_decimal!(
            df(3.0 * PI / 2.0, &[]),
            test_utils::cot_deriv(3.0 * PI / 2.0),
            15
        );
        assert_eq!(
            df(7.0 * PI / 4.0, &[]),
            test_utils::cot_deriv(7.0 * PI / 4.0)
        );
    }

    #[test]
    fn test_sderivative_inverse_sine() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.asin()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-0.5, &[]), test_utils::asin_deriv(-0.5));
        assert_eq!(df(0.0, &[]), test_utils::asin_deriv(0.0));
        assert_eq!(df(0.5, &[]), test_utils::asin_deriv(0.5));
    }

    #[test]
    fn test_sderivative_inverse_cosine() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.acos()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-0.5, &[]), test_utils::acos_deriv(-0.5));
        assert_eq!(df(0.0, &[]), test_utils::acos_deriv(0.0));
        assert_eq!(df(0.5, &[]), test_utils::acos_deriv(0.5));
    }

    #[test]
    fn test_sderivative_inverse_tangent() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.atan()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-1.5, &[]), test_utils::atan_deriv(-1.5));
        assert_eq!(df(-1.0, &[]), test_utils::atan_deriv(-1.0));
        assert_eq!(df(-0.5, &[]), test_utils::atan_deriv(-0.5));
        assert_eq!(df(0.0, &[]), test_utils::atan_deriv(0.0));
        assert_eq!(df(0.5, &[]), test_utils::atan_deriv(0.5));
        assert_eq!(df(1.0, &[]), test_utils::atan_deriv(1.0));
        assert_eq!(df(1.5, &[]), test_utils::atan_deriv(1.5));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_inverse_cosecant() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.acsc()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-1.5, &[]), test_utils::acsc_deriv(-1.5));
        assert_eq!(df(1.5, &[]), test_utils::acsc_deriv(1.5));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_inverse_secant() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.asec()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-1.5, &[]), test_utils::asec_deriv(-1.5));
        assert_eq!(df(1.5, &[]), test_utils::asec_deriv(1.5));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_inverse_cotangent() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.acot()
        }
        get_sderivative!(f, df);
        assert_equal_to_decimal!(df(-1.5, &[]), test_utils::acot_deriv(-1.5), 16);
        assert_eq!(df(-1.0, &[]), test_utils::acot_deriv(-1.0));
        assert_eq!(df(-0.5, &[]), test_utils::acot_deriv(-0.5));
        assert_eq!(df(0.5, &[]), test_utils::acot_deriv(0.5));
        assert_eq!(df(1.0, &[]), test_utils::acot_deriv(1.0));
        assert_equal_to_decimal!(df(1.5, &[]), test_utils::acot_deriv(1.5), 16);
    }

    #[test]
    fn test_sderivative_hyperbolic_sine() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.sinh()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-1.0, &[]), test_utils::sinh_deriv(-1.0));
        assert_eq!(df(0.0, &[]), test_utils::sinh_deriv(0.0));
        assert_eq!(df(1.0, &[]), test_utils::sinh_deriv(1.0));
    }

    #[test]
    fn test_sderivative_hyperbolic_cosine() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.cosh()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-1.0, &[]), test_utils::cosh_deriv(-1.0));
        assert_eq!(df(0.0, &[]), test_utils::cosh_deriv(0.0));
        assert_eq!(df(1.0, &[]), test_utils::cosh_deriv(1.0));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_hyperbolic_tangent() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.tanh()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-1.0, &[]), test_utils::tanh_deriv(-1.0));
        assert_eq!(df(0.0, &[]), test_utils::tanh_deriv(0.0));
        assert_eq!(df(1.0, &[]), test_utils::tanh_deriv(1.0));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_hyperbolic_cosecant() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.csch()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-1.0, &[]), test_utils::csch_deriv(-1.0));
        assert_eq!(df(1.0, &[]), test_utils::csch_deriv(1.0));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_hyperbolic_secant() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.sech()
        }
        get_sderivative!(f, df);
        assert_equal_to_decimal!(df(-1.0, &[]), test_utils::sech_deriv(-1.0), 16);
        assert_eq!(df(0.0, &[]), test_utils::sech_deriv(0.0));
        assert_equal_to_decimal!(df(1.0, &[]), test_utils::sech_deriv(1.0), 16);
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_hyperbolic_cotangent() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.coth()
        }
        get_sderivative!(f, df);
        assert_equal_to_decimal!(df(-1.0, &[]), test_utils::coth_deriv(-1.0), 16);
        assert_equal_to_decimal!(df(1.0, &[]), test_utils::coth_deriv(1.0), 16);
    }

    #[test]
    fn test_sderivative_inverse_hyperbolic_sine() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.asinh()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-1.5, &[]), test_utils::asinh_deriv(-1.5));
        assert_eq!(df(-1.0, &[]), test_utils::asinh_deriv(-1.0));
        assert_eq!(df(-0.5, &[]), test_utils::asinh_deriv(-0.5));
        assert_eq!(df(0.0, &[]), test_utils::asinh_deriv(0.0));
        assert_eq!(df(0.5, &[]), test_utils::asinh_deriv(0.5));
        assert_eq!(df(1.0, &[]), test_utils::asinh_deriv(1.0));
        assert_eq!(df(1.5, &[]), test_utils::asinh_deriv(1.5));
    }

    #[test]
    fn test_sderivative_inverse_hyperbolic_cosine() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.acosh()
        }
        get_sderivative!(f, df);
        assert_eq!(df(1.5, &[]), test_utils::acosh_deriv(1.5));
    }

    #[test]
    fn test_sderivative_inverse_hyperbolic_tangent() {
        fn f<S: Scalar>(x: S, _p: &[f64]) -> S {
            x.atanh()
        }
        get_sderivative!(f, df);
        assert_eq!(df(-0.5, &[]), test_utils::atanh_deriv(-0.5));
        assert_eq!(df(0.0, &[]), test_utils::atanh_deriv(0.0));
        assert_eq!(df(0.5, &[]), test_utils::atanh_deriv(0.5));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_inverse_hyperbolic_cosecant() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.acsch()
        }
        get_sderivative!(f, df);
        assert_equal_to_decimal!(df(-1.5, &[]), test_utils::acsch_deriv(-1.5), 16);
        assert_eq!(df(-1.0, &[]), test_utils::acsch_deriv(-1.0));
        assert_eq!(df(-0.5, &[]), test_utils::acsch_deriv(-0.5));
        assert_eq!(df(0.5, &[]), test_utils::acsch_deriv(0.5));
        assert_eq!(df(1.0, &[]), test_utils::acsch_deriv(1.0));
        assert_equal_to_decimal!(df(1.5, &[]), test_utils::acsch_deriv(1.5), 16);
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_inverse_hyperbolic_secant() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.asech()
        }
        get_sderivative!(f, df);
        assert_eq!(df(0.5, &[]), test_utils::asech_deriv(0.5));
    }

    #[test]
    #[cfg(feature = "trig")]
    fn test_sderivative_inverse_hyperbolic_cotangent() {
        fn f<S: Scalar + Trig>(x: S, _p: &[f64]) -> S {
            x.acoth()
        }
        get_sderivative!(f, df);
        assert_equal_to_decimal!(df(-1.5, &[]), test_utils::acoth_deriv(-1.5), 16);
        assert_equal_to_decimal!(df(1.5, &[]), test_utils::acoth_deriv(1.5), 16);
    }

    #[test]
    fn test_sderivative_with_runtime_parameters() {
        // Function to take the derivative of.
        fn f<S: Scalar>(x: S, p: &[f64]) -> S {
            let alpha = S::new(p[0]);
            let beta = S::new(p[1]);
            (alpha * x).exp() * (beta * x).cos()
        }

        // True derivative function.
        fn df<S: Scalar>(x: S, p: &[f64]) -> S {
            let alpha = S::new(p[0]);
            let beta = S::new(p[1]);
            let exp_term = (alpha * x).exp();
            let cos_term = (beta * x).cos();
            let sin_term = (beta * x).sin();
            exp_term * (alpha * cos_term - beta * sin_term)
        }

        // Evaluation point.
        let x0 = 0.3;

        // Parameter vector.
        let p = [0.5, 2.0];

        // Derivative function obtained via forward-mode automatic differentation.
        get_sderivative!(f, df_autodiff);

        // Evaluate the derivative using both functions.
        let df_eval_autodiff: f64 = df_autodiff(x0, &p);
        let df_eval: f64 = df(x0, &p);

        // Test autodiff derivative against true derivative.
        assert_equal_to_decimal!(df_eval_autodiff, df_eval, 16);
    }

    #[test]
    fn test_sderivative_custom_params() {
        struct Data {
            a: f64,
            b: f64,
            c: f64,
        }

        // Function to take the derivative of.
        #[allow(clippy::many_single_char_names)]
        fn f<S: Scalar>(x: S, p: &Data) -> S {
            let a = S::new(p.a);
            let b = S::new(p.b);
            let c = S::new(p.c);
            a * x.powi(2) + b * x + c
        }

        // Runtime parameter struct.
        let p = Data {
            a: 2.5,
            b: -1.3,
            c: 4.7,
        };

        // Derivative function obtained via forward-mode automatic differentiation.
        get_sderivative!(f, df, Data);

        // True derivative function.
        let df_true = |x: f64| 2.0 * p.a * x + p.b;

        // Evaluation point.
        let x0 = 1.0;

        // Evaluate the derivative using both functions.
        let df_eval: f64 = df(x0, &p);
        let df_eval_true: f64 = df_true(x0);

        // Test autodiff derivative against true derivative.
        assert_equal_to_decimal!(df_eval, df_eval_true, 15);
    }
}
