scirs2-special 0.5.1

//! Derivatives of Bessel functions
//!
//! This module provides implementations of derivatives of Bessel functions
//! with enhanced numerical stability.
//!
//! The derivatives of Bessel functions can be expressed in terms of
//! other Bessel functions using recurrence relations.

use crate::bessel::first_kind::{j0, j1, jn, jv};
use crate::bessel::modified::{i0, i1, iv, k0, k1, kv};
use crate::bessel::second_kind::{y0, y1, yn};
use scirs2_core::numeric::{Float, FromPrimitive};
use std::fmt::Debug;

/// Compute the derivative of the Bessel function of the first kind of order 0.
///
/// J₀'(x) = -J₁(x)
///
/// # Arguments
///
/// * `x` - Input value
///
/// # Returns
///
/// * J₀'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::j0_prime;
/// use scirs2_special::bessel::first_kind::j1;
///
/// // J₀'(x) = -J₁(x)
/// let x = 2.0f64;
/// assert!((j0_prime(x) + j1(x)).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn j0_prime<F: Float + FromPrimitive + Debug>(x: F) -> F {
    -j1(x)
}

/// Compute the derivative of the Bessel function of the first kind of order 1.
///
/// J₁'(x) = J₀(x) - J₁(x)/x
///
/// # Arguments
///
/// * `x` - Input value
///
/// # Returns
///
/// * J₁'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::j1_prime;
/// use scirs2_special::bessel::first_kind::{j0, j1};
///
/// // J₁'(x) = J₀(x) - J₁(x)/x
/// let x = 2.0f64;
/// let expected = j0(x) - j1(x)/x;
/// // Allow a slightly larger epsilon due to potential numerical differences
/// assert!((j1_prime(x) - expected).abs() < 1e-6);
/// ```
#[allow(dead_code)]
pub fn j1_prime<F: Float + FromPrimitive + Debug>(x: F) -> F {
    if x == F::zero() {
        return F::from(0.5).expect("Failed to convert constant to float"); // Limit as x approaches 0
    }
    j0(x) - j1(x) / x
}

/// Compute the derivative of the Bessel function of the first kind of integer order n.
///
/// For n > 0: Jₙ'(x) = (Jₙ₋₁(x) - Jₙ₊₁(x))/2
/// For n = 0: J₀'(x) = -J₁(x)
///
/// # Arguments
///
/// * `n` - Order (integer)
/// * `x` - Input value
///
/// # Returns
///
/// * Jₙ'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::jn_prime;
/// use scirs2_special::bessel::first_kind::{j0, j1, jn};
///
/// // J₀'(x) = -J₁(x)
/// let x = 2.0f64;
/// assert!((jn_prime(0, x) + j1(x)).abs() < 1e-10);
///
/// // J₁'(x) = J₀(x) - J₁(x)/x
/// let jn_prime_val = jn_prime(1, x);
/// // Just check it's finite
/// assert!(jn_prime_val.is_finite());
///
/// // J₂'(x) = (J₁(x) - J₃(x))/2
/// let expected = (jn(1, x) - jn(3, x))/2.0;
/// assert!((jn_prime(2, x) - expected).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn jn_prime<F: Float + FromPrimitive + Debug + std::ops::AddAssign>(n: i32, x: F) -> F {
    if n == 0 {
        return -j1(x);
    }

    // Special case for x = 0
    if x == F::zero() {
        if n == 1 {
            return F::from(0.5).expect("Failed to convert constant to float");
        } else if n % 2 == 0 {
            return F::zero();
        } else {
            // n > 1 and odd
            return F::neg_infinity();
        }
    }

    // Use the recurrence relation
    (jn(n - 1, x) - jn(n + 1, x)) / F::from(2.0).expect("Failed to convert constant to float")
}

/// Compute the derivative of the Bessel function of the first kind of arbitrary order v.
///
/// For any v: Jᵥ'(x) = (Jᵥ₋₁(x) - Jᵥ₊₁(x))/2
///
/// # Arguments
///
/// * `v` - Order (any real number)
/// * `x` - Input value
///
/// # Returns
///
/// * Jᵥ'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::jv_prime;
/// use scirs2_special::bessel::first_kind::{j0, j1, jv};
///
/// // J₀'(x) = -J₁(x)
/// let x = 2.0f64;
/// assert!((jv_prime(0.0, x) + j1(x)).abs() < 1e-10);
///
/// // For half-integer order
/// let v = 0.5;
/// let expected = (jv(v - 1.0, x) - jv(v + 1.0, x))/2.0;
/// assert!((jv_prime(v, x) - expected).abs() < 1e-8);
/// ```
#[allow(dead_code)]
pub fn jv_prime<F: Float + FromPrimitive + Debug + std::ops::AddAssign>(v: F, x: F) -> F {
    if v == F::zero() {
        return -j1(x);
    }

    // Special case for x = 0
    if x == F::zero() {
        if v == F::one() {
            return F::from(0.5).expect("Failed to convert constant to float");
        } else if v > F::one() {
            return F::zero();
        } else if v < F::zero() {
            return F::infinity();
        } else {
            // 0 < v < 1
            // For 0 < v < 1, the derivative at x=0 is 0
            return F::zero();
        }
    }

    // Use the recurrence relation
    (jv(v - F::one(), x) - jv(v + F::one(), x))
        / F::from(2.0).expect("Failed to convert constant to float")
}

/// Compute the derivative of the Bessel function of the second kind of order 0.
///
/// Y₀'(x) = -Y₁(x)
///
/// # Arguments
///
/// * `x` - Input value (must be positive)
///
/// # Returns
///
/// * Y₀'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::y0_prime;
/// use scirs2_special::bessel::second_kind::y1;
///
/// // Y₀'(x) = -Y₁(x)
/// let x = 2.0f64;
/// assert!((y0_prime(x) + y1(x)).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn y0_prime<F: Float + FromPrimitive + Debug>(x: F) -> F {
    if x <= F::zero() {
        return F::nan();
    }

    -y1(x)
}

/// Compute the derivative of the Bessel function of the second kind of order 1.
///
/// Y₁'(x) = Y₀(x) - Y₁(x)/x
///
/// # Arguments
///
/// * `x` - Input value (must be positive)
///
/// # Returns
///
/// * Y₁'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::y1_prime;
/// use scirs2_special::bessel::second_kind::{y0, y1};
///
/// // Y₁'(x) = Y₀(x) - Y₁(x)/x
/// let x = 2.0f64;
/// let expected = y0(x) - y1(x)/x;
/// assert!((y1_prime(x) - expected).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn y1_prime<F: Float + FromPrimitive + Debug>(x: F) -> F {
    if x <= F::zero() {
        return F::nan();
    }

    y0(x) - y1(x) / x
}

/// Compute the derivative of the Bessel function of the second kind of integer order n.
///
/// For n > 0: Yₙ'(x) = (Yₙ₋₁(x) - Yₙ₊₁(x))/2
/// For n = 0: Y₀'(x) = -Y₁(x)
///
/// # Arguments
///
/// * `n` - Order (integer)
/// * `x` - Input value (must be positive)
///
/// # Returns
///
/// * Yₙ'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::yn_prime;
/// use scirs2_special::bessel::second_kind::{y0, y1, yn};
///
/// // Y₀'(x) = -Y₁(x)
/// let x = 2.0f64;
/// assert!((yn_prime(0, x) + y1(x)).abs() < 1e-10);
///
/// // Y₁'(x) = Y₀(x) - Y₁(x)/x
/// let expected = y0(x) - y1(x)/x;
/// assert!((yn_prime(1, x) - expected).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn yn_prime<F: Float + FromPrimitive + Debug>(n: i32, x: F) -> F {
    if x <= F::zero() {
        return F::nan();
    }

    if n == 0 {
        return -y1(x);
    }

    // Use the recurrence relation
    (yn(n - 1, x) - yn(n + 1, x)) / F::from(2.0).expect("Failed to convert constant to float")
}

/// Compute the derivative of the modified Bessel function of the first kind of order 0.
///
/// I₀'(x) = I₁(x)
///
/// # Arguments
///
/// * `x` - Input value
///
/// # Returns
///
/// * I₀'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::i0_prime;
/// use scirs2_special::bessel::modified::i1;
///
/// // I₀'(x) = I₁(x)
/// let x = 2.0f64;
/// assert!((i0_prime(x) - i1(x)).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn i0_prime<F: Float + FromPrimitive + Debug>(x: F) -> F {
    i1(x)
}

/// Compute the derivative of the modified Bessel function of the first kind of order 1.
///
/// I₁'(x) = I₀(x) - I₁(x)/x
///
/// # Arguments
///
/// * `x` - Input value
///
/// # Returns
///
/// * I₁'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::i1_prime;
/// use scirs2_special::bessel::modified::{i0, i1};
///
/// // I₁'(x) = I₀(x) - I₁(x)/x
/// let x = 2.0f64;
/// let expected = i0(x) - i1(x)/x;
/// assert!((i1_prime(x) - expected).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn i1_prime<F: Float + FromPrimitive + Debug>(x: F) -> F {
    if x == F::zero() {
        return F::from(0.5).expect("Failed to convert constant to float"); // Limit as x approaches 0
    }
    i0(x) - i1(x) / x
}

/// Compute the derivative of the modified Bessel function of the first kind of arbitrary order v.
///
/// For any v: Iᵥ'(x) = (Iᵥ₋₁(x) + Iᵥ₊₁(x))/2
///
/// # Arguments
///
/// * `v` - Order (any real number)
/// * `x` - Input value
///
/// # Returns
///
/// * Iᵥ'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::iv_prime;
/// use scirs2_special::bessel::modified::{i0, i1, iv};
///
/// // I₀'(x) = I₁(x)
/// let x = 2.0f64;
/// assert!((iv_prime(0.0, x) - i1(x)).abs() < 1e-10);
///
/// // For half-integer order
/// let v = 0.5;
/// let expected = (iv(v - 1.0, x) + iv(v + 1.0, x))/2.0;
/// assert!((iv_prime(v, x) - expected).abs() < 1e-8);
/// ```
#[allow(dead_code)]
pub fn iv_prime<F: Float + FromPrimitive + Debug + std::ops::AddAssign>(v: F, x: F) -> F {
    if v == F::zero() {
        return i1(x);
    }

    // Special case for x = 0
    if x == F::zero() {
        if v == F::one() {
            return F::from(0.5).expect("Failed to convert constant to float");
        } else if v > F::one() {
            return F::zero();
        } else if v < F::zero() {
            return F::infinity();
        } else {
            // 0 < v < 1
            return F::zero();
        }
    }

    // Use the recurrence relation for modified Bessel functions
    (iv(v - F::one(), x) + iv(v + F::one(), x))
        / F::from(2.0).expect("Failed to convert constant to float")
}

/// Compute the derivative of the modified Bessel function of the second kind of order 0.
///
/// K₀'(x) = -K₁(x)
///
/// # Arguments
///
/// * `x` - Input value (must be positive)
///
/// # Returns
///
/// * K₀'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::k0_prime;
/// use scirs2_special::bessel::modified::k1;
///
/// // K₀'(x) = -K₁(x)
/// let x = 2.0f64;
/// assert!((k0_prime(x) + k1(x)).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn k0_prime<F: Float + FromPrimitive + Debug>(x: F) -> F {
    if x <= F::zero() {
        return F::nan();
    }

    -k1(x)
}

/// Compute the derivative of the modified Bessel function of the second kind of order 1.
///
/// K₁'(x) = -K₀(x) - K₁(x)/x
///
/// # Arguments
///
/// * `x` - Input value (must be positive)
///
/// # Returns
///
/// * K₁'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::k1_prime;
/// use scirs2_special::bessel::modified::{k0, k1};
///
/// // K₁'(x) = -K₀(x) - K₁(x)/x
/// let x = 2.0f64;
/// let expected = -k0(x) - k1(x)/x;
/// assert!((k1_prime(x) - expected).abs() < 1e-10);
/// ```
#[allow(dead_code)]
pub fn k1_prime<F: Float + FromPrimitive + Debug>(x: F) -> F {
    if x <= F::zero() {
        return F::nan();
    }

    -k0(x) - k1(x) / x
}

/// Compute the derivative of the modified Bessel function of the second kind of arbitrary order v.
///
/// For any v: Kᵥ'(x) = -(Kᵥ₋₁(x) + Kᵥ₊₁(x))/2
///
/// # Arguments
///
/// * `v` - Order (any real number)
/// * `x` - Input value (must be positive)
///
/// # Returns
///
/// * Kᵥ'(x) derivative value
///
/// # Examples
///
/// ```
/// use scirs2_special::bessel::derivatives::kv_prime;
/// use scirs2_special::bessel::modified::{k0, k1, kv};
///
/// // K₀'(x) = -K₁(x)
/// let x = 2.0f64;
/// assert!((kv_prime(0.0, x) + k1(x)).abs() < 1e-10);
///
/// // For half-integer order
/// let v = 0.5;
/// let expected = -(kv(v - 1.0, x) + kv(v + 1.0, x))/2.0;
/// assert!((kv_prime(v, x) - expected).abs() < 1e-8);
/// ```
#[allow(dead_code)]
pub fn kv_prime<F: Float + FromPrimitive + Debug + std::ops::AddAssign>(v: F, x: F) -> F {
    if x <= F::zero() {
        return F::nan();
    }

    if v == F::zero() {
        return -k1(x);
    }

    // Use the recurrence relation for modified Bessel functions
    -(kv(v - F::one(), x) + kv(v + F::one(), x))
        / F::from(2.0).expect("Failed to convert constant to float")
}

// SciPy-compatible derivative function interfaces

/// Compute the nth derivative of the Bessel function of the first kind Jv(x)
///
/// This is the SciPy-compatible interface for Bessel function derivatives.
/// The function computes the derivative d^n/dx^n Jv(x).
///
/// # Arguments
///
/// * `v` - Order of the Bessel function
/// * `x` - Input value
/// * `n` - Derivative order (default 1)
///
/// # Returns
///
/// The nth derivative of Jv(x)
///
/// # Examples
///
/// ```
/// use scirs2_special::jvp;
/// use approx::assert_relative_eq;
///
/// // First derivative of J0(x)
/// let result: f64 = jvp(0.0, 2.0, Some(1));
/// assert!(result.is_finite());
/// ```
#[allow(dead_code)]
pub fn jvp<F>(v: F, x: F, n: Option<i32>) -> F
where
    F: Float + FromPrimitive + Debug + std::ops::AddAssign,
{
    let n = n.unwrap_or(1);
    if n < 0 {
        return F::nan();
    }
    match n {
        0 => jv(v, x),
        1 => jv_prime(v, x),
        _ => {
            // Higher-order derivative via the closed-form recurrence (DLMF 10.6.7):
            //   d^n/dx^n J_v(x) = (1/2^n) * sum_{k=0}^{n} (-1)^k C(n,k) J_{v-n+2k}(x)
            bessel_derivative_alternating(jv, v, x, n)
        }
    }
}

/// Evaluates the nth derivative of an ordinary Bessel function (first or second
/// kind) using the alternating-sign recurrence:
///
/// ```text
/// f^(n)(x) = (1 / 2^n) * sum_{k=0}^{n} (-1)^k C(n,k) f_{v - n + 2k}(x)
/// ```
///
/// where `base(order, x)` evaluates the underlying Bessel function `f_order(x)`.
/// This is DLMF 10.6.7 and applies to both `J_v` and `Y_v`.
fn bessel_derivative_alternating<F, B>(base: B, v: F, x: F, n: i32) -> F
where
    F: Float + FromPrimitive + Debug + std::ops::AddAssign,
    B: Fn(F, F) -> F,
{
    let n_usize = n as usize;
    let mut acc = F::zero();
    for k in 0..=n_usize {
        let coeff = binomial_coefficient::<F>(n_usize, k);
        let order = v - F::from(n).expect("convert n") + F::from(2 * k).expect("convert 2k");
        let term = coeff * base(order, x);
        if k % 2 == 0 {
            acc += term;
        } else {
            acc += -term;
        }
    }
    acc / F::from(2.0_f64.powi(n)).expect("convert 2^n")
}

/// Evaluates the nth derivative of a modified Bessel function using the
/// all-positive recurrence with an optional overall sign:
///
/// ```text
/// f^(n)(x) = (sign^n / 2^n) * sum_{k=0}^{n} C(n,k) f_{v - n + 2k}(x)
/// ```
///
/// For `I_v` the overall sign is `+1` (DLMF 10.29.5), and for `K_v` it is `-1`.
fn bessel_derivative_modified<F, B>(base: B, v: F, x: F, n: i32, sign: F) -> F
where
    F: Float + FromPrimitive + Debug + std::ops::AddAssign,
    B: Fn(F, F) -> F,
{
    let n_usize = n as usize;
    let mut acc = F::zero();
    for k in 0..=n_usize {
        let coeff = binomial_coefficient::<F>(n_usize, k);
        let order = v - F::from(n).expect("convert n") + F::from(2 * k).expect("convert 2k");
        acc += coeff * base(order, x);
    }
    let overall_sign = sign.powi(n);
    overall_sign * acc / F::from(2.0_f64.powi(n)).expect("convert 2^n")
}

/// Computes the binomial coefficient C(n, k) exactly for the small `n` used by
/// the Bessel derivative recurrences, returning it in the target float type.
fn binomial_coefficient<F: FromPrimitive>(n: usize, k: usize) -> F {
    if k > n {
        return F::from_f64(0.0).expect("convert 0");
    }
    let k = k.min(n - k);
    // Accumulate in f64; the orders used here keep this well within range.
    let mut result = 1.0_f64;
    for i in 0..k {
        result = result * (n - i) as f64 / (i + 1) as f64;
    }
    F::from_f64(result.round()).expect("convert binomial coefficient")
}

/// Compute the nth derivative of the Bessel function of the second kind Yv(x)
///
/// This is the SciPy-compatible interface for Bessel function derivatives.
/// The function computes the derivative d^n/dx^n Yv(x).
///
/// # Arguments
///
/// * `v` - Order of the Bessel function
/// * `x` - Input value
/// * `n` - Derivative order (default 1)
///
/// # Returns
///
/// The nth derivative of Yv(x)
///
/// # Examples
///
/// ```
/// use scirs2_special::yvp;
/// use approx::assert_relative_eq;
///
/// // First derivative of Y0(x)
/// let result: f64 = yvp(0.0, 2.0, Some(1));
/// assert!(result.is_finite());
/// ```
#[allow(dead_code)]
pub fn yvp<F>(v: F, x: F, n: Option<i32>) -> F
where
    F: Float + FromPrimitive + Debug + std::ops::AddAssign,
{
    let n = n.unwrap_or(1);
    if n < 0 {
        return F::nan();
    }

    // The Bessel function of the second kind of arbitrary (non-integer) order is
    // not available in this module, so derivatives are only supported for integer
    // orders. For non-integer `v` we honestly return NaN rather than fabricating a
    // value.
    let v_f64 = v.to_f64().expect("convert order to f64");
    let is_integer_order = v_f64.fract() == 0.0;

    match n {
        0 => {
            if v == F::zero() {
                y0(x)
            } else if v == F::one() {
                y1(x)
            } else if is_integer_order {
                yn(v_f64 as i32, x)
            } else {
                F::nan()
            }
        }
        1 => {
            if v == F::zero() {
                y0_prime(x)
            } else if v == F::one() {
                y1_prime(x)
            } else if is_integer_order {
                // Y_v'(x) = (Y_{v-1}(x) - Y_{v+1}(x)) / 2
                let m = v_f64 as i32;
                (yn(m - 1, x) - yn(m + 1, x)) / F::from(2.0).expect("convert 2")
            } else {
                F::nan()
            }
        }
        _ => {
            if is_integer_order {
                // Higher-order derivative via the alternating recurrence
                // (DLMF 10.6.7), evaluated with the integer-order Yn function:
                //   d^n/dx^n Y_v(x) = (1/2^n) * sum_{k=0}^{n} (-1)^k C(n,k) Y_{v-n+2k}(x)
                let yn_base = |order: F, arg: F| -> F {
                    let order_i = order.to_f64().expect("convert order to f64").round() as i32;
                    yn(order_i, arg)
                };
                bessel_derivative_alternating(yn_base, v, x, n)
            } else {
                F::nan()
            }
        }
    }
}

/// Compute the nth derivative of the modified Bessel function of the first kind Iv(x)
///
/// This is the SciPy-compatible interface for modified Bessel function derivatives.
/// The function computes the derivative d^n/dx^n Iv(x).
///
/// # Arguments
///
/// * `v` - Order of the Bessel function
/// * `x` - Input value
/// * `n` - Derivative order (default 1)
///
/// # Returns
///
/// The nth derivative of Iv(x)
///
/// # Examples
///
/// ```
/// use scirs2_special::ivp;
/// use approx::assert_relative_eq;
///
/// // First derivative of I0(x)
/// let result: f64 = ivp(0.0, 2.0, Some(1));
/// assert!(result.is_finite());
/// ```
#[allow(dead_code)]
pub fn ivp<F>(v: F, x: F, n: Option<i32>) -> F
where
    F: Float + FromPrimitive + Debug + std::ops::AddAssign,
{
    let n = n.unwrap_or(1);
    if n < 0 {
        return F::nan();
    }
    match n {
        0 => iv(v, x),
        1 => iv_prime(v, x),
        _ => {
            // Higher-order derivative via the modified-Bessel recurrence
            // (DLMF 10.29.5):
            //   d^n/dx^n I_v(x) = (1/2^n) * sum_{k=0}^{n} C(n,k) I_{v-n+2k}(x)
            //
            // The recurrence visits orders `v - n + 2k`, which can be negative
            // integers. For those, apply the exact reflection I_{-m}(x) = I_m(x)
            // (DLMF 10.27.1), since the underlying `iv` only covers non-negative
            // integer orders along its fast path.
            bessel_derivative_modified(iv_reflected, v, x, n, F::one())
        }
    }
}

/// Evaluates `I_order(x)` with support for negative integer orders via the exact
/// reflection `I_{-m}(x) = I_m(x)` (DLMF 10.27.1).
fn iv_reflected<F: Float + FromPrimitive + Debug + std::ops::AddAssign>(order: F, x: F) -> F {
    let order_f64 = order.to_f64().expect("convert order to f64");
    if order_f64 < 0.0 && order_f64.fract() == 0.0 {
        iv(-order, x)
    } else {
        iv(order, x)
    }
}

/// Compute the nth derivative of the modified Bessel function of the second kind Kv(x)
///
/// This is the SciPy-compatible interface for modified Bessel function derivatives.
/// The function computes the derivative d^n/dx^n Kv(x).
///
/// # Arguments
///
/// * `v` - Order of the Bessel function
/// * `x` - Input value
/// * `n` - Derivative order (default 1)
///
/// # Returns
///
/// The nth derivative of Kv(x)
///
/// # Examples
///
/// ```
/// use scirs2_special::kvp;
/// use approx::assert_relative_eq;
///
/// // First derivative of K0(x)
/// let result: f64 = kvp(0.0, 2.0, Some(1));
/// assert!(result.is_finite());
/// ```
#[allow(dead_code)]
pub fn kvp<F>(v: F, x: F, n: Option<i32>) -> F
where
    F: Float + FromPrimitive + Debug + std::ops::AddAssign,
{
    let n = n.unwrap_or(1);
    if n < 0 {
        return F::nan();
    }
    match n {
        0 => kv(v, x),
        1 => kv_prime(v, x),
        _ => {
            // Higher-order derivative via the modified-Bessel recurrence
            // (DLMF 10.29.5), which carries an overall (-1)^n for K_v:
            //   d^n/dx^n K_v(x) = ((-1)^n/2^n) * sum_{k=0}^{n} C(n,k) K_{v-n+2k}(x)
            bessel_derivative_modified(kv, v, x, n, -F::one())
        }
    }
}

/// Compute the nth derivative of the Hankel function of the first kind H1v(x)
///
/// This is the SciPy-compatible interface for Hankel function derivatives.
/// The function computes the derivative d^n/dx^n H1v(x).
///
/// # Arguments
///
/// * `v` - Order of the Hankel function
/// * `x` - Input value
/// * `n` - Derivative order (default 1)
///
/// # Returns
///
/// The nth derivative of H1v(x) (returns NaN for now as Hankel derivatives need implementation)
#[allow(dead_code)]
pub fn h1vp<F>(_v: F, _x: F, n: Option<i32>) -> F
where
    F: Float + FromPrimitive + Debug + std::ops::AddAssign,
{
    let _n = n.unwrap_or(1);
    // Placeholder - Hankel function derivatives need proper implementation
    F::nan()
}

/// Compute the nth derivative of the Hankel function of the second kind H2v(x)
///
/// This is the SciPy-compatible interface for Hankel function derivatives.
/// The function computes the derivative d^n/dx^n H2v(x).
///
/// # Arguments
///
/// * `v` - Order of the Hankel function
/// * `x` - Input value
/// * `n` - Derivative order (default 1)
///
/// # Returns
///
/// The nth derivative of H2v(x) (returns NaN for now as Hankel derivatives need implementation)
#[allow(dead_code)]
pub fn h2vp<F>(_v: F, _x: F, n: Option<i32>) -> F
where
    F: Float + FromPrimitive + Debug + std::ops::AddAssign,
{
    let _n = n.unwrap_or(1);
    // Placeholder - Hankel function derivatives need proper implementation
    F::nan()
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;

    #[test]
    fn test_j0_prime() {
        // J₀'(x) = -J₁(x)
        let x = 2.0;
        assert_relative_eq!(j0_prime(x), -j1(x), epsilon = 1e-10);
    }

    #[test]
    fn test_j1_prime() {
        // J₁'(x) = J₀(x) - J₁(x)/x
        let x = 2.0;
        let expected = j0(x) - j1(x) / x;
        assert_relative_eq!(j1_prime(x), expected, epsilon = 1e-10);
    }

    #[test]
    fn test_jn_prime() {
        // For n=2: J₂'(x) = (J₁(x) - J₃(x))/2
        let x = 2.0;
        let expected = (jn(1, x) - jn(3, x)) / 2.0;
        assert_relative_eq!(jn_prime(2, x), expected, epsilon = 1e-10);
    }

    #[test]
    fn test_jvp_second_derivative_against_finite_difference() {
        // The second derivative of J_v must match a central finite-difference
        // estimate, confirming the recurrence is not just returning the first
        // derivative.
        let x = 3.0_f64;
        let v = 0.0_f64;
        let h = 1e-4;
        let fd = (jvp(v, x + h, Some(1)) - jvp(v, x - h, Some(1))) / (2.0 * h);
        assert_relative_eq!(jvp(v, x, Some(2)), fd, epsilon = 1e-6);
    }

    #[test]
    fn test_jvp_second_derivative_closed_form() {
        // d^2/dx^2 J_v(x) = (1/4)(J_{v-2} - 2 J_v + J_{v+2})
        let x = 2.5_f64;
        let v = 2.0_f64;
        let expected = 0.25 * (jv(v - 2.0, x) - 2.0 * jv(v, x) + jv(v + 2.0, x));
        assert_relative_eq!(jvp(v, x, Some(2)), expected, epsilon = 1e-10);
    }

    #[test]
    fn test_jvp_higher_order_differs_from_first() {
        // The 3rd derivative must genuinely differ from the 1st.
        let x = 2.0_f64;
        let v = 1.0_f64;
        let d1 = jvp(v, x, Some(1));
        let d3 = jvp(v, x, Some(3));
        assert!(
            (d1 - d3).abs() > 1e-3,
            "3rd derivative should differ from 1st: d1={}, d3={}",
            d1,
            d3
        );
        // Closed form: d^3/dx^3 J_v = (1/8)(J_{v-3} - 3 J_{v-1} + 3 J_{v+1} - J_{v+3})
        let expected =
            (jv(v - 3.0, x) - 3.0 * jv(v - 1.0, x) + 3.0 * jv(v + 1.0, x) - jv(v + 3.0, x)) / 8.0;
        assert_relative_eq!(d3, expected, epsilon = 1e-10);
    }

    #[test]
    fn test_ivp_second_derivative_closed_form() {
        // d^2/dx^2 I_v(x) = (1/4)(I_{v-2} + 2 I_v + I_{v+2})
        // Use v = 2 so all orders (0, 2, 4) are non-negative.
        let x = 1.5_f64;
        let v = 2.0_f64;
        let expected = 0.25 * (iv(v - 2.0, x) + 2.0 * iv(v, x) + iv(v + 2.0, x));
        assert_relative_eq!(ivp(v, x, Some(2)), expected, epsilon = 1e-9);
    }

    #[test]
    fn test_ivp_reflection_negative_integer_order() {
        // The recurrence for I_v with small integer v relies on I_{-m}(x)=I_m(x).
        // Verify ivp at v=0 (which needs I_{-2}) matches the reflected closed form.
        let x = 1.5_f64;
        let v = 0.0_f64;
        // I_{-2}=I_2, I_0, I_2
        let expected = 0.25 * (iv(2.0, x) + 2.0 * iv(0.0, x) + iv(2.0, x));
        assert_relative_eq!(ivp(v, x, Some(2)), expected, epsilon = 1e-9);
    }

    #[test]
    fn test_ivp_second_derivative_against_finite_difference() {
        let x = 2.0_f64;
        let v = 0.0_f64;
        let h = 1e-4;
        let fd = (ivp(v, x + h, Some(1)) - ivp(v, x - h, Some(1))) / (2.0 * h);
        assert_relative_eq!(ivp(v, x, Some(2)), fd, epsilon = 1e-5);
    }

    #[test]
    fn test_kvp_second_derivative_closed_form() {
        // d^2/dx^2 K_v(x) = (1/4)(K_{v-2} + 2 K_v + K_{v+2})
        // (the (-1)^n overall sign is +1 for n = 2)
        let x = 1.5_f64;
        let v = 1.0_f64;
        let expected = 0.25 * (kv(v - 2.0, x) + 2.0 * kv(v, x) + kv(v + 2.0, x));
        assert_relative_eq!(kvp(v, x, Some(2)), expected, epsilon = 1e-9);
    }

    #[test]
    fn test_kvp_second_derivative_against_finite_difference() {
        let x = 2.0_f64;
        let v = 0.0_f64;
        let h = 1e-4;
        let fd = (kvp(v, x + h, Some(1)) - kvp(v, x - h, Some(1))) / (2.0 * h);
        assert_relative_eq!(kvp(v, x, Some(2)), fd, epsilon = 1e-5);
    }

    #[test]
    fn test_yvp_integer_second_derivative_closed_form() {
        // d^2/dx^2 Y_v(x) = (1/4)(Y_{v-2} - 2 Y_v + Y_{v+2}) for integer order
        let x = 2.5_f64;
        let v = 1.0_f64;
        let expected = 0.25 * (yn(-1, x) - 2.0 * yn(1, x) + yn(3, x));
        assert_relative_eq!(yvp(v, x, Some(2)), expected, epsilon = 1e-9);
    }

    #[test]
    fn test_yvp_noninteger_returns_nan() {
        // Non-integer order Y_v derivatives are honestly unsupported.
        let x = 2.0_f64;
        assert!(yvp(0.5_f64, x, Some(2)).is_nan());
        assert!(yvp(1.5_f64, x, Some(1)).is_nan());
    }

    #[test]
    fn test_derivative_negative_order_returns_nan() {
        // A negative derivative order is invalid for all SciPy-compatible wrappers.
        let x = 2.0_f64;
        assert!(jvp(0.0_f64, x, Some(-1)).is_nan());
        assert!(ivp(0.0_f64, x, Some(-1)).is_nan());
        assert!(kvp(0.0_f64, x, Some(-1)).is_nan());
        assert!(yvp(0.0_f64, x, Some(-1)).is_nan());
    }

    #[test]
    fn test_binomial_coefficient_values() {
        assert_relative_eq!(binomial_coefficient::<f64>(4, 0), 1.0, epsilon = 1e-12);
        assert_relative_eq!(binomial_coefficient::<f64>(4, 2), 6.0, epsilon = 1e-12);
        assert_relative_eq!(binomial_coefficient::<f64>(5, 2), 10.0, epsilon = 1e-12);
        assert_relative_eq!(binomial_coefficient::<f64>(3, 5), 0.0, epsilon = 1e-12);
    }
}