optionstratlib 0.17.3

/******************************************************************************
   Author: Joaquín Béjar García
   Email: jb@taunais.com
   Date: 11/8/24
******************************************************************************/

use crate::Options;
use crate::error::decimal::DecimalError;
use crate::error::greeks::{DeltaNeutralityErrorKind, GreeksError, InputErrorKind, MathErrorKind};
use crate::model::decimal::{d_div, d_mul, d_sub, f64_to_decimal};
use crate::strategies::DELTA_THRESHOLD;
use core::f64;
use num_traits::ToPrimitive;
use positive::Positive;
use rust_decimal::{Decimal, MathematicalOps};
use rust_decimal_macros::dec;
use statrs::distribution::{ContinuousCDF, Normal};

/// Calculates the `d1` parameter used in the Black-Scholes options pricing model.
///
/// The `d1` value is an intermediary result used to determine option greeks and prices.
/// It is computed using the formula:
///
/// ```math
/// d1 = (ln(S / K) + (r + σ² / 2) * T) / (σ * sqrt(T))
/// ```
///
/// Where:
/// - `S`: Underlying price
/// - `K`: Strike price
/// - `r`: Risk-free rate
/// - `T`: Time to expiration (in years)
/// - `σ`: Implied volatility
///
/// # Parameters
///
/// - `underlying_price`: The current price of the underlying asset. Must be positive.
/// - `strike_price`: The strike price of the option. Must be greater than zero.
/// - `risk_free_rate`: The annual risk-free interest rate, expressed as a decimal.
/// - `expiration_date`: The time to expiration of the option, in years. Must be greater than zero.
/// - `implied_volatility`: The implied volatility of the option, expressed as a decimal. Must be greater than zero.
///
/// # Returns
///
/// - `Ok(Decimal)`: The computed `d1` value.
/// - `Err(GreeksError)`: Returns an error if input validation fails. Possible errors include:
///   - Invalid underlying price (must be greater than zero).
///   - Invalid strike price (must be greater than zero).
///   - Invalid implied volatility (must be greater than zero).
///   - Invalid expiration time (must be greater than zero).
///
/// # Errors
///
/// Returns a `GreeksError::InputError` in the following cases:
/// - **InvalidPrice**: Triggered when `underlying_price` is zero or less.
/// - **InvalidStrike**: Triggered when `strike_price` is zero or less.
/// - **InvalidVolatility**: Triggered when `implied_volatility` is zero.
/// - **InvalidTime**: Triggered when `expiration_date` is zero or less.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::d1;
/// use positive::{pos_or_panic, Positive};
///
/// let underlying_price = Positive::HUNDRED;
/// let strike_price = pos_or_panic!(95.0);
/// let risk_free_rate = dec!(0.05);
/// let expiration_date = pos_or_panic!(0.5); // 6 months
/// let implied_volatility = pos_or_panic!(0.2);
///
/// match d1(
///     underlying_price,
///     strike_price,
///     risk_free_rate,
///     expiration_date,
///     implied_volatility,
/// ) {
///     Ok(result) => info!("d1: {}", result),
///     Err(e) => error!("Error: {:?}", e),
/// }
/// ```
#[inline]
pub fn d1(
    underlying_price: Positive,
    strike_price: Positive,
    risk_free_rate: Decimal,
    expiration_date: Positive,
    implied_volatility: Positive,
) -> Result<Decimal, GreeksError> {
    if underlying_price == Positive::ZERO {
        return Err(GreeksError::InputError(InputErrorKind::InvalidPrice {
            value: underlying_price.to_f64(),
            reason: "Underlying price price cannot be zero".to_string(),
        }));
    }

    if strike_price == Positive::ZERO {
        return Err(GreeksError::InputError(InputErrorKind::InvalidStrike {
            value: strike_price.to_string(),
            reason: "Strike price cannot be zero".to_string(),
        }));
    }

    if implied_volatility == Decimal::ZERO {
        return Err(GreeksError::InputError(InputErrorKind::InvalidVolatility {
            value: implied_volatility.to_f64(),
            reason: "Implied volatility cannot be zero".to_string(),
        }));
    }

    if expiration_date == Decimal::ZERO {
        return Err(GreeksError::InputError(InputErrorKind::InvalidTime {
            value: expiration_date,
            reason: "Expiration date cannot be zero".to_string(),
        }));
    }

    // d1 = (ln(S / K) + (r + σ² / 2) * T) / (σ * sqrt(T))
    let underlying_price: Decimal = underlying_price.to_dec();
    let implied_volatility_squared = implied_volatility.powd(Decimal::TWO);
    let ln_price_ratio = match strike_price {
        value if value == Positive::INFINITY => Decimal::MIN,
        _ => (underlying_price / strike_price).ln(),
    };

    let rate_vol_term = risk_free_rate + implied_volatility_squared / Decimal::TWO;
    let numerator = ln_price_ratio + rate_vol_term * expiration_date;
    let denominator = implied_volatility * expiration_date.sqrt();

    match numerator.checked_div(denominator.into()) {
        Some(result) => Ok(result),
        None => Err(GreeksError::MathError(MathErrorKind::Overflow)),
    }
}

/// Calculates the `d2` parameter used in the Black-Scholes options pricing model.
///
/// The `d2` value is an intermediary result derived from the `d1` value and is used
/// to determine option greeks and prices. It is computed using the formula:
///
/// ```math
/// d2 = d1 - σ * sqrt(T)
/// ```
///
/// Where:
/// - `d1`: The `d1` value calculated using the `d1` function.
/// - `σ`: Implied volatility.
/// - `T`: Time to expiration (in years).
///
/// # Parameters
///
/// - `underlying_price`: The current price of the underlying asset. Must be positive.
/// - `strike_price`: The strike price of the option. Must be greater than zero.
/// - `risk_free_rate`: The annual risk-free interest rate, expressed as a decimal.
/// - `expiration_date`: The time to expiration of the option, in years. Must be greater than zero.
/// - `implied_volatility`: The implied volatility of the option, expressed as a decimal. Must be greater than zero.
///
/// # Returns
///
/// - `Ok(Decimal)`: The computed `d2` value.
/// - `Err(GreeksError)`: Returns an error if input validation fails or if the `d1` computation fails.
///
/// # Errors
///
/// Returns a `GreeksError::InputError` in the following cases:
/// - **InvalidVolatility**: Triggered when `implied_volatility` is zero.
/// - **InvalidTime**: Triggered when `expiration_date` is zero.
///
/// # Notes
///
/// This function depends on the `d1` function to compute the `d1` value. Any errors from
/// the `d1` function will propagate to this function.
///
/// # Example
///
/// ```rust
/// # fn run() -> Result<(), optionstratlib::error::Error> {
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::d2;
/// use positive::{pos_or_panic, Positive};
/// let underlying_price = Positive::new(100.0)?;
/// let strike_price = Positive::new(95.0)?;
/// let risk_free_rate = dec!(0.05);
/// let expiration_date = pos_or_panic!(0.5); // 6 months
/// let implied_volatility = pos_or_panic!(0.2);
///
/// match d2(
///     underlying_price,
///     strike_price,
///     risk_free_rate,
///     expiration_date,
///     implied_volatility,
/// ) {
///     Ok(result) => info!("d2: {}", result),
///     Err(e) => error!("Error: {:?}", e),
/// }
/// # Ok(())
/// # }
/// ```
#[inline]
pub fn d2(
    underlying_price: Positive,
    strike_price: Positive,
    risk_free_rate: Decimal,
    expiration_date: Positive,
    implied_volatility: Positive,
) -> Result<Decimal, GreeksError> {
    //    function d1 is already checking for validity of implied_volatility and
    //    expiration_date with error propagation. Can we comment out the two if blocks below?
    //    if implied_volatility == Decimal::ZERO {
    //        return Err(GreeksError::InputError(InputErrorKind::InvalidVolatility {
    //            value: implied_volatility.to_f64(),
    //            reason: "Implied volatility cannot be zero".to_string(),
    //        }));
    //    }
    //
    //    if expiration_date == Decimal::ZERO {
    //        return Err(GreeksError::InputError(InputErrorKind::InvalidTime {
    //            value: expiration_date,
    //            reason: "Expiration date cannot be zero".to_string(),
    //        }));
    //    }

    let d1_value = d1(
        underlying_price,
        strike_price,
        risk_free_rate,
        expiration_date,
        implied_volatility,
    )?;

    Ok(d1_value - implied_volatility * expiration_date.sqrt())
}

/// Computes the probability density function (PDF) of the standard normal distribution
/// for a given input `x`.
///
/// The PDF of the standard normal distribution is defined as:
///
/// ```math
/// N(x) = \frac{1}{\sqrt{2 \pi}} \cdot e^{-\frac{x^2}{2}}
/// ```
///
/// Where:
/// - \(x\): The input value for which the PDF is computed.
///
/// # Parameters
///
/// - `x: Decimal`
///   The input value for which the standard normal PDF is calculated.
///
/// # Returns
///
/// - `Ok(Decimal)`: The computed PDF value as a `Decimal`.
/// - `Err(GreeksError)`: Returns an error if the computation fails.
///
/// # Calculation Details
///
/// - The normalisation factor \(1 / \sqrt{2 \pi}\) is held as a precomputed
///   `Decimal` constant (`dec!(0.3989422804014326779399461)`), avoiding a
///   runtime fallible `sqrt()` on `Decimal`.
/// - The exponent is computed as \(-\frac{x^2}{2}\).
/// - The PDF value is the product of the normalisation factor and the
///   exponential term.
/// - For very-negative exponents (`pre_pdf < -11.7`) `exp` underflows to
///   `0` in `Decimal` arithmetic, so we short-circuit to `Decimal::ZERO`.
///
/// # Errors
///
/// - `GreeksError`: This function will return an error if any part of the calculation fails,
///   though this is unlikely as the operations are well-defined for all finite inputs.
///
/// # Example
///
/// ```rust
/// use rust_decimal::Decimal;
/// use tracing::{error, info};
/// use optionstratlib::greeks::n;
///
/// let x = Decimal::new(100, 2); // 1.00
///
/// match n(x) {
///     Ok(result) => info!("N(x): {}", result),
///     Err(e) => error!("Error calculating N(x): {:?}", e),
/// }
/// ```
///
/// # Notes
///
/// The function uses the `Decimal` type for precision and error handling. The result is returned
#[inline]
pub fn n(x: Decimal) -> Result<Decimal, GreeksError> {
    // 1 / sqrt(2π) — standard normal PDF normalisation constant.
    // Precomputed so we avoid the runtime fallible sqrt on Decimal.
    const NORM_FACTOR: Decimal = dec!(0.3989422804014326779399461);
    let pre_pdf = -x.powd(Decimal::TWO) / Decimal::TWO;

    // avoid Exp underflowed
    if pre_pdf < dec!(-11.7) {
        return Ok(Decimal::ZERO);
    }

    let pdf = pre_pdf.exp();
    Ok(NORM_FACTOR * pdf) // N(x) = [1 / sqrt(2 * PI)] * e^(-x^2 / 2)
}

/// Calculate the derivative of the function `n` at a given point `x`.
///
/// This function represents the negative product of `x` and the result of the
/// function `n` evaluated at `x`.
///
/// # Arguments
///
/// * `x` - A floating point number representing the input to the function `n`.
///
/// # Returns
///
/// This function returns a value of type `T` which is the derivative of the
/// function `n` at `x`.
///
/// # Type Parameters
///
/// * `T`: A type that implements the `Float` trait, which allows for floating point operations.
///
/// # Examples
///
/// Although no example usage is provided, it typically involves calling `n_prime`
/// with a specific floating point number as the argument to obtain the derivative
/// of the function `n` at that point.
///
/// # Panics
///
/// This function does not explicitly handle panics. Ensure that the function `n`
/// and the floating point operations are well-defined for the input value `x`.
///
/// # Safety
///
/// This code assumes that the function `n` is defined and behaves correctly for
/// the provided input type `T`. Improper or undefined behavior of `n` may lead
/// to unexpected results or runtime errors.
#[allow(dead_code)]
#[inline]
pub(crate) fn n_prime(x: Decimal) -> Result<Decimal, GreeksError> {
    Ok(-x * n(x)?) // -x * n(x)
}

/// Computes the cumulative distribution function (CDF) of the standard normal distribution
/// for a given input `x`.
///
/// The function uses the standard normal distribution (mean = 0, standard deviation = 1)
/// to calculate the probability that a normally distributed random variable is less than or
/// equal to `x`. This is commonly referred to as `N(x)` in financial and statistical contexts.
///
/// # Parameters
///
/// - `x: Decimal`
///   The input value for which the CDF is computed. Must be convertible to `f64`.
///
/// # Returns
///
/// - `Ok(Decimal)`: The CDF value corresponding to the input `x`.
/// - `Err(DecimalError)`: Returns an error if the conversion from `Decimal` to `f64` fails.
///
/// # Errors
///
/// Returns a `DecimalError::ConversionError` if:
/// - The input `x` cannot be converted to an `f64`.
///
/// # Notes
///
/// This function uses the [`statrs`](https://docs.rs/statrs/latest/statrs/) crate to model the
/// standard normal distribution and compute the CDF. The result is returned as a `Decimal`
/// for precision.
///
/// # Example
///
/// ```rust
/// use rust_decimal::Decimal;
/// use tracing::{error, info};
/// use optionstratlib::greeks::big_n;
///
/// let x = Decimal::new(100, 2);
///
/// match big_n(x) {
///     Ok(result) => info!("N(x): {}", result),
///     Err(e) => error!("Error: {:?}", e),
/// }
/// ```
#[inline]
pub fn big_n(x: Decimal) -> Result<Decimal, DecimalError> {
    let Some(x_f64) = x.to_f64() else {
        return Err(DecimalError::ConversionError {
            from_type: "Decimal".to_string(),
            to_type: "f64".to_string(),
            reason: "Conversion failed".to_string(),
        });
    };

    // Guard the `Decimal` → `f64` boundary: if `x` is outside the
    // representable `f64` range the conversion returns `±∞` rather
    // than `None`, and feeding that into `statrs::cdf` yields `NaN`
    // which would later collapse silently in `f64_to_decimal`.
    if !x_f64.is_finite() {
        return Err(DecimalError::invalid_value(
            x_f64,
            "big_n: Decimal -> f64 produced a non-finite value",
        ));
    }

    const MEAN: f64 = 0.0;
    const STD_DEV: f64 = 1.0;

    // Normal::new(0.0, 1.0) is infallible by construction (parameters
    // hardcoded), but surface the error type via map_err to satisfy
    // the panic-free rule.
    let normal_distribution =
        Normal::new(MEAN, STD_DEV).map_err(|e| DecimalError::ConversionError {
            from_type: "(f64, f64)".to_string(),
            to_type: "Normal".to_string(),
            reason: format!("invalid Normal parameters: {e}"),
        })?;
    let cdf = normal_distribution.cdf(x_f64);
    if !cdf.is_finite() {
        return Err(DecimalError::invalid_value(
            cdf,
            "big_n: CDF produced a non-finite value",
        ));
    }
    f64_to_decimal(cdf)
}

/// Calculates the d1 and d2 values used in financial option pricing models such as the Black-Scholes model.
///
/// # Arguments
///
/// * `option` - A reference to an `Options` struct containing the underlying price,
///   the risk-free rate, and the implied volatility of the option.
/// * `time_to_expiry` - The time remaining until option expiry in years.
///
/// # Returns
///
/// * A tuple containing two `f64` values:
///     - `d1_value`: The calculated d1 value.
///     - `d2_value`: The calculated d2 value.
///
#[inline]
pub(crate) fn calculate_d_values(option: &Options) -> Result<(Decimal, Decimal), GreeksError> {
    let b = option.risk_free_rate - option.dividend_yield.to_dec();
    let d1_value = d1(
        option.underlying_price,
        option.strike_price,
        b,
        option.expiration_date.get_years()?,
        option.implied_volatility,
    );
    let d2_value = d2(
        option.underlying_price,
        option.strike_price,
        b,
        option.expiration_date.get_years()?,
        option.implied_volatility,
    );
    Ok((d1_value?, d2_value?))
}

/// Calculates d1 and d2 for the Black-76 model (options on futures/forwards).
///
/// For Black-76, the cost of carry `b = 0` because the forward price F already
/// incorporates all carry. This function specializes `calculate_d_values` for
/// that case.
///
/// # Arguments
/// * `option` - The option data; `underlying_price` holds the forward price F.
///
/// # Returns
/// * `Ok((d1, d2))` - The Black-76 d1 and d2 parameters.
/// * `Err(GreeksError)` - If validation or computation fails.
pub fn calculate_d_values_black_76(option: &Options) -> Result<(Decimal, Decimal), GreeksError> {
    // Black-76: cost of carry b = 0 (forward pricing, no carry term in d1/d2)
    let b = Decimal::ZERO;
    let years = option.expiration_date.get_years()?;
    let d1_value = d1(
        option.underlying_price,
        option.strike_price,
        b,
        years,
        option.implied_volatility,
    );
    let d2_value = d2(
        option.underlying_price,
        option.strike_price,
        b,
        years,
        option.implied_volatility,
    );
    Ok((d1_value?, d2_value?))
}

/// Calculates the optimal position sizes for two positions to achieve delta neutrality
/// while maintaining a specified total position size.
///
/// # Arguments
/// * `delta1` - Delta of the first position (e.g., short call delta)
/// * `delta2` - Delta of the second position (e.g., short put delta)
/// * `total_size` - Desired total position size (sum of both positions)
///
/// # Returns
/// * `Ok((size1, size2))` - Tuple containing the calculated sizes for each position
/// * `Err(String)` - Error message if calculation is not possible
///
/// # Example
/// ```rust
/// # fn run() -> Result<(), optionstratlib::error::Error> {
/// use rust_decimal_macros::dec;
/// use optionstratlib::greeks::calculate_delta_neutral_sizes;
/// use positive::pos_or_panic;
/// let (call_size, put_size) = calculate_delta_neutral_sizes(
///     dec!(-0.30),  // Short call delta
///     dec!(0.20),   // Short put delta
///     pos_or_panic!(7.0)     // Total desired position size
/// )?;
/// # Ok(())
/// # }
/// ```
///
/// # Errors
///
/// Returns [`GreeksError::DeltaNeutrality`] with
/// [`DeltaNeutralityErrorKind::ZeroDelta`] when both deltas are zero,
/// [`DeltaNeutralityErrorKind::EqualDeltas`] when `delta1 == delta2`,
/// [`DeltaNeutralityErrorKind::SameSignDeltas`] when the two deltas
/// share the same sign (no neutral mix exists), or
/// [`DeltaNeutralityErrorKind::NegativePositionSize`] when the derived
/// sizes fall below zero for the provided `total_size`.
pub fn calculate_delta_neutral_sizes(
    delta1: Decimal,
    delta2: Decimal,
    total_size: Positive,
) -> Result<(Positive, Positive), GreeksError> {
    // The equation we want to solve is:
    // delta1 * size1 + delta2 * size2 = Decimal::ZERO
    // size1 + size2 = total_size

    if delta1.is_zero() || delta2.is_zero() {
        return Err(DeltaNeutralityErrorKind::ZeroDelta.into());
    }

    // Validate inputs
    if delta1 == delta2 {
        return Err(DeltaNeutralityErrorKind::EqualDeltas.into());
    }

    if delta1.is_sign_positive() == delta2.is_sign_positive() {
        return Err(DeltaNeutralityErrorKind::SameSignDeltas.into());
    }

    // Delta-neutral sizing: size1 = -total_size · delta2 / (delta1 - delta2).
    // The multiplication and subtraction are both monetary (signed position
    // quantities on a potentially large total_size) and the division is the
    // only site where the outcome is projected onto the `Positive` invariant,
    // so overflow here has to surface explicitly rather than silently wrap.
    let weighted = d_mul(
        -total_size.to_dec(),
        delta2,
        "greeks::delta_neutral::size1::weighted",
    )?;
    let spread = d_sub(delta1, delta2, "greeks::delta_neutral::size1::spread")?;
    let size1 = Positive::new_decimal(d_div(weighted, spread, "greeks::delta_neutral::size1")?)?;
    let size2 = total_size - size1;

    // Validate results
    if size1 < Decimal::ZERO || size2 < Decimal::ZERO {
        return Err(DeltaNeutralityErrorKind::NegativePositionSize.into());
    }

    // Verify the solution
    let total_delta: Decimal = size1.to_dec() * delta1 + size2.to_dec() * delta2;
    if total_delta.abs() > DELTA_THRESHOLD {
        // Allow small numerical errors
        return Err(DeltaNeutralityErrorKind::NotAchievable.into());
    }
    let total_size_check = size1 + size2;
    if (total_size_check.to_dec() - total_size.to_dec()).abs() > DELTA_THRESHOLD {
        return Err(DeltaNeutralityErrorKind::SizeMismatch {
            calculated: total_size_check,
            expected: total_size,
        }
        .into());
    }

    Ok((size1, size2))
}

#[cfg(test)]
mod tests_calculate_delta_neutral_sizes {
    use super::*;
    use crate::assert_decimal_eq;
    use positive::{assert_pos_relative_eq, pos_or_panic};
    use rust_decimal_macros::dec;

    #[test]
    fn test_valid_delta_neutral_calculation() {
        let result = calculate_delta_neutral_sizes(
            dec!(-0.30),        // Short call delta
            dec!(0.20),         // Short put delta
            pos_or_panic!(7.0), // Total size
        )
        .unwrap();

        let (size1, size2) = result;

        // Check total size constraint
        assert_pos_relative_eq!(size1 + size2, pos_or_panic!(7.0), pos_or_panic!(0.0001));

        // Check delta neutrality (mixed-sign deltas — compute in Decimal).
        let total_delta = size1.to_dec() * dec!(-0.30) + size2.to_dec() * dec!(0.20);
        assert_decimal_eq!(total_delta, Decimal::ZERO, dec!(0.0001));
    }

    #[test]
    fn test_equal_deltas() {
        let result = calculate_delta_neutral_sizes(dec!(0.25), dec!(0.25), pos_or_panic!(10.0));
        assert!(result.is_err());
    }

    #[test]
    fn test_impossible_neutrality() {
        let result = calculate_delta_neutral_sizes(dec!(-0.95), dec!(-0.90), pos_or_panic!(10.0));
        assert!(result.is_err());
    }
}

#[cfg(test)]
mod tests_calculate_d_values {
    use super::*;
    use crate::model::types::{OptionStyle, OptionType, Side};

    use approx::assert_relative_eq;
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    #[test]
    fn test_calculate_d_values() {
        let option = Options {
            option_type: OptionType::European,
            side: Side::Long,
            underlying_symbol: "".to_string(),
            strike_price: pos_or_panic!(110.0),
            underlying_price: Positive::HUNDRED,
            risk_free_rate: dec!(0.05),
            implied_volatility: pos_or_panic!(10.12),
            expiration_date: Default::default(),
            quantity: Positive::ONE,
            option_style: OptionStyle::Call,
            dividend_yield: Positive::ZERO,
            exotic_params: None,
        };
        let (d1_value, d2_value) = calculate_d_values(&option).unwrap();

        assert_relative_eq!(
            d1_value.to_f64().unwrap(),
            5.055522709505501,
            epsilon = 0.001
        );
        assert_relative_eq!(
            d2_value.to_f64().unwrap(),
            -5.064477290494499,
            epsilon = 0.001
        );
    }
}

#[cfg(test)]
mod tests_src_greeks_utils {
    use super::*;

    use approx::assert_relative_eq;
    use num_traits::FloatConst;
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;
    use statrs::distribution::ContinuousCDF;
    use statrs::distribution::Normal;

    #[test]
    fn test_d1_zero_sigma() {
        let s = Positive::HUNDRED;
        let k = Positive::HUNDRED;
        let r = dec!(0.05);
        let t = Positive::ONE;
        let sigma = Positive::ZERO;
        let _ = d1(s, k, r, t, sigma).is_err();
    }

    #[test]
    fn test_d1_zero_t() {
        let s = Positive::HUNDRED;
        let k = Positive::HUNDRED;
        let r = dec!(0.05);
        let t = Positive::ZERO;
        let sigma = pos_or_panic!(0.01);
        let _ = d1(s, k, r, t, sigma).is_err();
    }

    #[test]
    fn test_d2_bis_i() {
        let s = Positive::HUNDRED;
        let k = pos_or_panic!(110.0);
        let r = dec!(0.05);
        let t = Positive::TWO;
        let sigma = pos_or_panic!(0.2);
        let computed_d2 = d2(s, k, r, t, sigma).unwrap().to_f64().unwrap();
        let computed_d1 = d1(s, k, r, t, sigma).unwrap().to_f64().unwrap();
        assert_relative_eq!(computed_d1, 0.15800237455184707, epsilon = 0.001);
        assert_relative_eq!(computed_d2, -0.12484033792277195, epsilon = 0.001);
    }

    #[test]
    fn test_d2_bis_ii() {
        let s = Positive::HUNDRED;
        let k = pos_or_panic!(95.0);
        let r = dec!(0.15);
        let t = Positive::ONE;
        let sigma = pos_or_panic!(0.2);
        let computed_d2 = d2(s, k, r, t, sigma).unwrap().to_f64().unwrap();
        let computed_d1 = d1(s, k, r, t, sigma).unwrap().to_f64().unwrap();
        assert_relative_eq!(computed_d1, 1.1064664719377526, epsilon = 0.001);
        assert_relative_eq!(computed_d2, 0.9064664719377528, epsilon = 0.001);
    }

    #[test]
    fn test_d2_zero_sigma() {
        let s = Positive::HUNDRED;
        let k = Positive::HUNDRED;
        let r = Decimal::ZERO;
        let t = Positive::ONE;
        let sigma = Positive::ZERO;
        let _ = d2(s, k, r, t, sigma).is_err();
    }

    #[test]
    fn test_d2_zero_t() {
        let s = Positive::HUNDRED;
        let k = Positive::HUNDRED;
        let r = dec!(0.02);
        let t = Positive::ZERO;
        let sigma = pos_or_panic!(0.01);
        let _ = d2(s, k, r, t, sigma).is_err();
    }

    #[test]
    fn test_n() {
        let x = Decimal::ZERO;
        let expected_n = 1.0 / (2.0 * f64::PI()).sqrt();
        let computed_n = n(x).unwrap().to_f64().unwrap();
        assert_relative_eq!(computed_n, expected_n, epsilon = 1e-8);

        let x = Decimal::ONE;
        let expected_n = 1.0 / (2.0 * f64::PI()).sqrt() * (-0.5f64).exp();
        let computed_n = n(x).unwrap().to_f64().unwrap();
        assert_relative_eq!(computed_n, expected_n, epsilon = 1e-8);
    }

    #[test]
    fn test_big_n() {
        let x = Decimal::ZERO;
        let normal_distribution = Normal::new(0.0, 1.0).unwrap();
        let expected_big_n = normal_distribution.cdf(x.to_f64().unwrap());
        let computed_big_n = big_n(x).unwrap().to_f64().unwrap();
        assert!(
            (computed_big_n - expected_big_n).abs() < 1e-10,
            "big_n function failed"
        );

        let x = Decimal::ONE;
        let expected_big_n = normal_distribution.cdf(1.0);
        let computed_big_n = big_n(x).unwrap().to_f64().unwrap();
        assert!(
            (computed_big_n - expected_big_n).abs() < 1e-10,
            "big_n function failed"
        );
    }
}

#[cfg(test)]
mod calculate_d1_values {
    use super::*;
    use positive::pos_or_panic;

    use rust_decimal_macros::dec;

    #[test]
    fn test_d1_zero_volatility() {
        // Case where volatility (sigma) is zero
        let underlying_price = Positive::HUNDRED;
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ONE;
        let implied_volatility = Positive::ZERO;

        // When volatility is zero, d1 should handle the case correctly
        assert!(
            d1(
                underlying_price,
                strike_price,
                risk_free_rate,
                expiration_date,
                implied_volatility,
            )
            .is_err()
        );
    }

    #[test]
    fn test_d1_zero_time_to_expiry() {
        // Case where time to expiry is zero
        let underlying_price = Positive::HUNDRED;
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ZERO;
        let implied_volatility = pos_or_panic!(0.2);

        // When time to expiry is zero, d1 should handle the case correctly
        assert!(
            d1(
                underlying_price,
                strike_price,
                risk_free_rate,
                expiration_date,
                implied_volatility,
            )
            .is_err()
        );
    }

    #[test]
    fn test_d1_high_volatility() {
        // Case with extremely high volatility
        let underlying_price = Positive::HUNDRED;
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ONE;
        let implied_volatility = Positive::HUNDRED; // Very high volatility

        // High volatility should result in a small or large value for d1
        let calculated_d1 = d1(
            underlying_price,
            strike_price,
            risk_free_rate,
            expiration_date,
            implied_volatility,
        )
        .unwrap()
        .to_f64()
        .unwrap();

        // Assert the result should be finite and non-infinite
        assert!(
            calculated_d1.is_finite(),
            "d1 should not be infinite for high volatility"
        );
    }

    #[test]
    fn test_d1_high_underlying_price() {
        // Case with extremely high underlying price
        let underlying_price = Positive::INFINITY; // Very high stock price
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ONE;
        let implied_volatility = pos_or_panic!(0.2);

        // Very high underlying price should result in a large d1 value
        assert!(
            d1(
                underlying_price,
                strike_price,
                risk_free_rate,
                expiration_date,
                implied_volatility,
            )
            .is_ok()
        );
    }

    #[test]
    fn test_d1_low_underlying_price() {
        // Case with extremely low underlying price (near zero)
        let underlying_price = pos_or_panic!(0.01); // Very low stock price
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ONE;
        let implied_volatility = pos_or_panic!(0.2);

        // Very low underlying price should result in a small or negative d1 value
        let calculated_d1 = d1(
            underlying_price,
            strike_price,
            risk_free_rate,
            expiration_date,
            implied_volatility,
        )
        .unwrap()
        .to_f64()
        .unwrap();

        // Assert the result should be finite and not infinite
        assert!(
            calculated_d1.is_finite(),
            "d1 should not be infinite for low underlying price"
        );
    }

    #[test]
    fn test_d1_zero_strike_price() {
        // Case where strike price is zero
        let underlying_price = Positive::HUNDRED;
        let strike_price = Positive::ZERO;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ONE;
        let implied_volatility = pos_or_panic!(0.2);

        // Since strike price is zero, the function should call handle_zero and return positive infinity
        assert!(
            d1(
                underlying_price,
                strike_price,
                risk_free_rate,
                expiration_date,
                implied_volatility,
            )
            .is_err()
        );
    }

    #[test]
    fn test_d1_infinite_risk_free_rate() {
        // Case where risk-free rate is very high (infinite-like)
        let underlying_price = Positive::HUNDRED;
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = Decimal::MAX; // Very high risk-free rate
        let expiration_date = Positive::ONE;
        let implied_volatility = pos_or_panic!(0.2);

        // High risk-free rate should result in a large d1 value, potentially infinite
        assert!(
            d1(
                underlying_price,
                strike_price,
                risk_free_rate,
                expiration_date,
                implied_volatility,
            )
            .is_err()
        );
    }
}

#[cfg(test)]
mod calculate_d1_values_bis {
    use super::*;
    use crate::error::greeks::{GreeksError, InputErrorKind};

    use approx::assert_relative_eq;
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    // Helper function to convert Decimal to f64 for testing
    fn decimal_to_f64_test(d: Decimal) -> f64 {
        d.to_f64().unwrap()
    }

    #[test]
    fn test_d1_basic_calculation() {
        let result = d1(
            Positive::HUNDRED,
            pos_or_panic!(90.0),
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        );

        assert!(result.is_ok());
        let d1_value = decimal_to_f64_test(result.unwrap());
        assert_relative_eq!(d1_value, 0.8768025782891316, epsilon = 0.0001);
    }

    #[test]
    fn test_d1_in_the_money() {
        let result = d1(
            pos_or_panic!(110.0),
            pos_or_panic!(90.0),
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        );

        assert!(result.is_ok());
        let d1_value = decimal_to_f64_test(result.unwrap());
        assert_relative_eq!(d1_value, 1.3533534773107558, epsilon = 0.0001);
    }

    #[test]
    fn test_d1_out_of_the_money() {
        let result = d1(
            pos_or_panic!(90.0),
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        );

        assert!(result.is_ok());
        let d1_value = decimal_to_f64_test(result.unwrap());
        assert_relative_eq!(d1_value, -0.1768025782891315, epsilon = 0.0001);
    }

    #[test]
    fn test_d1_zero_strike_error() {
        let result = d1(
            Positive::HUNDRED,
            Positive::ZERO,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        );

        assert!(matches!(
            result,
            Err(GreeksError::InputError(
                InputErrorKind::InvalidStrike { .. }
            ))
        ));
    }

    #[test]
    fn test_d1_zero_volatility_error() {
        let result = d1(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            Positive::ZERO,
        );

        assert!(matches!(
            result,
            Err(GreeksError::InputError(
                InputErrorKind::InvalidVolatility { .. }
            ))
        ));
    }

    #[test]
    fn test_d1_zero_time_error() {
        let result = d1(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ZERO,
            pos_or_panic!(0.2),
        );

        assert!(matches!(
            result,
            Err(GreeksError::InputError(InputErrorKind::InvalidTime { .. }))
        ));
    }

    #[test]
    fn test_d1_short_expiry() {
        let result = d1(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            pos_or_panic!(0.0833), // approximately one month
            pos_or_panic!(0.05),
        );

        assert!(result.is_ok());
        let d1_value = decimal_to_f64_test(result.unwrap());
        assert_relative_eq!(d1_value, 0.29583282863806715, epsilon = 0.0001);
    }

    #[test]
    fn test_d1_high_volatility() {
        let result = d1(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.5), // 50% volatility
        );

        assert!(result.is_ok());
        let d1_value = decimal_to_f64_test(result.unwrap());
        assert_relative_eq!(d1_value, 0.35, epsilon = 0.0001);
    }

    #[test]
    fn test_d1_zero_interest_rate() {
        let result = d1(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.0),
            Positive::ONE,
            pos_or_panic!(0.5),
        );

        assert!(result.is_ok());
        let d1_value = decimal_to_f64_test(result.unwrap());
        assert_relative_eq!(d1_value, 0.25, epsilon = 0.0001);
    }

    #[test]
    fn test_d1_negative_interest_rate() {
        let result = d1(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(-0.02), // negative interest rate
            Positive::ONE,
            pos_or_panic!(0.5),
        );

        assert!(result.is_ok());
        let d1_value = decimal_to_f64_test(result.unwrap());
        assert_relative_eq!(d1_value, 0.21, epsilon = 0.0001);
    }

    #[test]
    fn test_d1_negative_interest_rate_bis() {
        let result = d1(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(-0.02), // negative interest rate
            Positive::ONE,
            pos_or_panic!(0.5),
        );

        assert!(result.is_ok());
        let d1_value = decimal_to_f64_test(result.unwrap());
        assert_relative_eq!(d1_value, 0.21, epsilon = 0.0001);
    }
}

#[cfg(test)]
mod calculate_d2_values {
    use super::*;
    use positive::pos_or_panic;

    use rust_decimal_macros::dec;

    #[test]
    fn test_d2_zero_volatility() {
        // Case where volatility (implied_volatility) is zero
        let underlying_price = Positive::HUNDRED;
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ONE;
        let implied_volatility = Positive::ZERO;

        // When volatility is zero, d2 should handle the case correctly using handle_zero
        assert!(
            d2(
                underlying_price,
                strike_price,
                risk_free_rate,
                expiration_date,
                implied_volatility,
            )
            .is_err()
        );
    }

    #[test]
    fn test_d2_zero_time_to_expiry() {
        // Case where time to expiration is zero
        let underlying_price = Positive::HUNDRED;
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ZERO;
        let implied_volatility = pos_or_panic!(0.2);

        // When time to expiration is zero, handle_zero should be called
        assert!(
            d2(
                underlying_price,
                strike_price,
                risk_free_rate,
                expiration_date,
                implied_volatility,
            )
            .is_err()
        );
    }

    #[test]
    fn test_d2_high_volatility() {
        // Case with extremely high volatility
        let underlying_price = Positive::HUNDRED;
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ONE;
        let implied_volatility = Positive::HUNDRED; // Very high volatility

        // High volatility should result in a significant negative shift in d2
        let calculated_d2 = d2(
            underlying_price,
            strike_price,
            risk_free_rate,
            expiration_date,
            implied_volatility,
        )
        .unwrap()
        .to_f64()
        .unwrap();

        // d2 should be finite and not infinite
        assert!(
            calculated_d2.is_finite(),
            "d2 should not be infinite for high volatility"
        );
    }

    #[test]
    fn test_d2_high_underlying_price() {
        // Case with extremely high underlying price
        let underlying_price = Positive::INFINITY;
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ONE;
        let implied_volatility = pos_or_panic!(0.2);

        // Very high underlying price should result in a large d2 value
        assert!(
            d2(
                underlying_price,
                strike_price,
                risk_free_rate,
                expiration_date,
                implied_volatility,
            )
            .is_ok()
        );
    }

    #[test]
    fn test_d2_low_underlying_price() {
        // Case with extremely low underlying price (near zero)
        let underlying_price = pos_or_panic!(0.01);
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ONE;
        let implied_volatility = pos_or_panic!(0.2);

        // Very low underlying price should result in a small or negative d2 value
        let calculated_d2 = d2(
            underlying_price,
            strike_price,
            risk_free_rate,
            expiration_date,
            implied_volatility,
        )
        .unwrap()
        .to_f64()
        .unwrap();

        // Assert the result should be finite and not infinite
        assert!(
            calculated_d2.is_finite(),
            "d2 should not be infinite for low underlying price"
        );
    }

    #[test]
    fn test_d2_zero_strike_price() {
        // Case where strike price is zero
        let underlying_price = Positive::HUNDRED;
        let strike_price = Positive::ZERO;
        let risk_free_rate = dec!(0.05);
        let expiration_date = Positive::ONE;
        let implied_volatility = pos_or_panic!(0.2);

        // Since strike price is zero, the function should call handle_zero and return positive infinity
        assert!(
            d2(
                underlying_price,
                strike_price,
                risk_free_rate,
                expiration_date,
                implied_volatility,
            )
            .is_err()
        );
    }

    #[test]
    fn test_d2_infinite_risk_free_rate() {
        // Case where risk-free rate is very high (infinite-like)
        let underlying_price = Positive::HUNDRED;
        let strike_price = Positive::HUNDRED;
        let risk_free_rate = Decimal::MAX; // Very high risk-free rate
        let expiration_date = Positive::ONE;
        let implied_volatility = pos_or_panic!(0.2);

        // High risk-free rate should result in a large d2 value, potentially infinite
        assert!(
            d2(
                underlying_price,
                strike_price,
                risk_free_rate,
                expiration_date,
                implied_volatility,
            )
            .is_err()
        );
    }
}

#[cfg(test)]
mod calculate_d2_values_bis {
    use super::*;
    use crate::assert_decimal_eq;
    use approx::assert_relative_eq;
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    const EPSILON: Decimal = dec!(0.0001);
    // Normal test cases
    #[test]
    fn test_d2_atm_option() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_relative_eq!(result.to_f64().unwrap(), 0.15, epsilon = 0.0001);
    }

    #[test]
    fn test_d2_itm_call() {
        let result = d2(
            pos_or_panic!(110.0),
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_decimal_eq!(result, dec!(0.6265508990216243), EPSILON);
    }

    #[test]
    fn test_d2_otm_call() {
        let result = d2(
            pos_or_panic!(90.0),
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            -0.3768025782891315,
            epsilon = 0.0001
        );
    }

    // Time to expiration variations
    #[test]
    fn test_d2_short_expiry() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            pos_or_panic!(0.0833), // 1 month
            pos_or_panic!(0.5),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            -0.04329260906898544,
            epsilon = 0.0001
        );
    }

    #[test]
    fn test_d2_long_expiry() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::TWO,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            0.21213203435596426,
            epsilon = 0.0001
        );
    }

    // Volatility variations
    #[test]
    fn test_d2_low_volatility() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.1),
        )
        .unwrap();
        assert_relative_eq!(result.to_f64().unwrap(), 0.45, epsilon = 0.0001);
    }

    #[test]
    fn test_d2_high_volatility() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.5),
        )
        .unwrap();
        assert_relative_eq!(result.to_f64().unwrap(), -0.15, epsilon = 0.0001);
    }

    // Interest rate variations
    #[test]
    fn test_d2_zero_interest() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            Decimal::ZERO,
            Positive::ONE,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_relative_eq!(result.to_f64().unwrap(), -0.1, epsilon = 0.0001);
    }

    #[test]
    fn test_d2_high_interest() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.10),
            Positive::ONE,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_relative_eq!(result.to_f64().unwrap(), 0.4, epsilon = 0.0001);
    }

    // Extreme price differences
    #[test]
    fn test_d2_deep_itm() {
        let result = d2(
            pos_or_panic!(200.0),
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            3.6157359027997265,
            epsilon = 0.0001
        );
    }

    #[test]
    fn test_d2_deep_otm() {
        let result = d2(
            pos_or_panic!(50.0),
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            -3.3157359027997266,
            epsilon = 0.0001
        );
    }

    // Very small values
    #[test]
    fn test_d2_small_price() {
        let result = d2(
            pos_or_panic!(0.01),
            pos_or_panic!(0.01),
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_relative_eq!(result.to_f64().unwrap(), 0.15, epsilon = 0.0001);
    }

    #[test]
    fn test_d2_small_time() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            pos_or_panic!(0.001),
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            0.004743416490252569,
            epsilon = 0.0001
        );
    }

    #[test]
    fn test_d2_small_volatility() {
        let result = d2(
            pos_or_panic!(200.0),
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.01),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            74.30971805599454,
            epsilon = 0.0001
        );
    }

    // Error cases
    #[test]
    fn test_d2_zero_volatility() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            Positive::ZERO,
        );
        assert!(matches!(
            result,
            Err(GreeksError::InputError(
                InputErrorKind::InvalidVolatility { .. }
            ))
        ));
    }

    #[test]
    fn test_d2_zero_time() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ZERO,
            pos_or_panic!(0.2),
        );
        assert!(matches!(
            result,
            Err(GreeksError::InputError(InputErrorKind::InvalidTime { .. }))
        ));
    }

    // Negative interest rate
    #[test]
    fn test_d2_negative_interest() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            -dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_decimal_eq!(result, dec!(-0.35), EPSILON);
    }

    // Combined extreme cases
    #[test]
    fn test_d2_combined_extremes_high() {
        let result = d2(
            pos_or_panic!(1000.0),
            Positive::HUNDRED,
            dec!(0.15),
            pos_or_panic!(5.0),
            pos_or_panic!(0.8),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            0.812019752759385,
            epsilon = 0.0001
        );
    }

    #[test]
    fn test_d2_combined_extremes_low() {
        let result = d2(
            pos_or_panic!(10.0),
            Positive::HUNDRED,
            dec!(0.01),
            pos_or_panic!(0.1),
            pos_or_panic!(0.05),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            -145.57292814518308,
            epsilon = 0.0001
        );
    }

    // Edge cases with very large numbers
    #[test]
    fn test_d2_large_price_ratio() {
        let result = d2(
            pos_or_panic!(1_000_000.0),
            Positive::ONE,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            69.22755278982137,
            epsilon = 0.0001
        );
    }

    // Special case: ATM LEAPS (Long-term equity anticipation securities)
    #[test]
    fn test_d2_leaps() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            pos_or_panic!(2.5), // 2.5 years
            pos_or_panic!(0.15),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            0.40846086443841567,
            epsilon = 0.0001
        );
    }

    // Near-zero but valid cases
    #[test]
    fn test_d2_near_zero_valid_values() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.0001),
            pos_or_panic!(0.01),
            pos_or_panic!(0.001),
        )
        .unwrap();
        assert!(result.to_f64().unwrap().abs() < 1.0);
    }

    // Test with maximum realistic market values
    #[test]
    fn test_d2_max_realistic_values() {
        let result = d2(
            pos_or_panic!(10000.0),
            pos_or_panic!(5000.0),
            dec!(0.20),
            pos_or_panic!(3.0),
            pos_or_panic!(1.5),
        )
        .unwrap();
        assert_relative_eq!(
            result.to_f64().unwrap(),
            -0.8013055238112647,
            epsilon = 0.0001
        );
    }
}

#[cfg(test)]
mod calculate_n_values {
    use super::*;
    use approx::assert_relative_eq;
    use rust_decimal_macros::dec;
    use std::f64::consts::PI;

    #[test]
    fn test_n_zero() {
        // Case where x = 0.0
        let x = Decimal::ZERO;

        // The PDF of the standard normal distribution at x = 0 is 1/sqrt(2*pi)
        let expected_n = 1.0f64 / (2.0 * PI).sqrt();

        // Compute n(x)
        let calculated_n = n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n, expected_n, epsilon = 1e-8);
    }

    #[test]
    fn test_n_one() {
        // Case where x = 0.0
        let x = Decimal::ONE;

        // The PDF of the standard normal distribution at x = 1 is 0.24197072535043143
        let expected_n = 0.24197072535043143;

        // Compute n(x)
        let calculated_n = n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n, expected_n, epsilon = 1e-8);
    }

    #[test]
    fn test_n_two() {
        // Case where x = 0.0
        let x = Decimal::TWO;

        // The PDF of the standard normal distribution at x = 2 is 0.05399096672219953
        let expected_n = 0.05399096672219953;

        // Compute n(x)
        let calculated_n = n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n, expected_n, epsilon = 1e-8);
    }

    #[test]
    fn test_n_positive_small_value() {
        // Case where x is a small positive value
        let x = dec!(0.5);

        // Expected result for n(0.5), can be precomputed
        let expected_n = 1.0f64 / (2.0 * PI).sqrt() * (-0.5f64 * 0.5f64 / 2.0f64).exp();

        // Compute n(x)
        let calculated_n = n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n, expected_n, epsilon = 1e-8);
    }

    #[test]
    fn test_n_negative_small_value() {
        // Case where x is a small negative value
        let x = dec!(-0.5);

        // Expected result for n(-0.5), which should be the same as n(0.5) due to symmetry
        let expected_n = 1.0f64 / (2.0 * PI).sqrt() * (-0.5f64 * 0.5f64 / 2.0f64).exp();

        // Compute n(x)
        let calculated_n = n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n, expected_n, epsilon = 1e-8);
    }

    #[test]
    fn test_n_large_positive_value() {
        // Case where x is a large positive value
        let x = dec!(5.0);

        // Compute n(x)
        let calculated_n = n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n, 0.0, epsilon = 1e-8);
    }

    #[test]
    fn test_n_large_negative_value() {
        // Case where x is a large negative value
        let x = dec!(-5.0);

        // Compute n(x)
        let calculated_n = n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n, 0.0, epsilon = 1e-8);
    }

    #[test]
    fn test_n_extreme_positive_value() {
        // Case where x is a very large positive value
        let x = dec!(100.0);

        // Expected result for n(100.0), should be extremely close to 0
        let expected_n = 1.0f64 / (2.0 * PI).sqrt() * (-100.0f64 * 100.0f64 / 2.0f64).exp();

        // Compute n(x)
        let calculated_n = n(x).unwrap().to_f64().unwrap();

        // Assert that n(x) is effectively 0 for such a large input
        assert_relative_eq!(calculated_n, expected_n, epsilon = 1e-100);
    }

    #[test]
    fn test_n_extreme_negative_value() {
        // Case where x is a very large negative value
        let x = dec!(-100.0);

        // Expected result for n(-100.0), should be extremely close to 0
        let expected_n = 1.0f64 / (2.0 * PI).sqrt() * (-100.0f64 * 100.0f64 / 2.0f64).exp();

        // Compute n(x)
        let calculated_n = n(x).unwrap().to_f64().unwrap();

        // Assert that n(x) is effectively 0 for such a large negative input
        assert_relative_eq!(calculated_n, expected_n, epsilon = 1e-100);
    }
}

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

    #[test]
    fn test_n_prime_zero() {
        // Case where x = 0.0
        let x = dec!(0.0);

        // The derivative of the PDF at x = 0 should be 0 because -x * n(x) = 0 * n(0) = 0
        let expected_n_prime = 0.0f64;

        // Compute n_prime(x)
        let calculated_n_prime = n_prime(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n_prime, expected_n_prime, epsilon = 1e-8);
    }

    #[test]
    fn test_n_prime_one() {
        // Case where x = 0.0
        let x = Decimal::ONE;

        // The derivative of the PDF at x = 1 should be 0.24197072535043143f64 (n(1) = 0.24197072535043143)
        let expected_n_prime = -0.24197072535043143f64;

        // Compute n_prime(x)
        let calculated_n_prime = n_prime(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n_prime, expected_n_prime, epsilon = 1e-8);
    }

    #[test]
    fn test_n_prime_two() {
        // Case where x = 0.0
        let x = Decimal::TWO;

        // The derivative of the PDF at x = 2 should be -0.10798193344439906 (n(2) = 0.05399096672219953)
        let expected_n_prime = -0.10798193344439906;

        // Compute n_prime(x)
        let calculated_n_prime = n_prime(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n_prime, expected_n_prime, epsilon = 1e-8);
    }

    #[test]
    fn test_n_prime_positive_small_value() {
        // Case where x is a small positive value
        let x = dec!(0.5);

        // Expected result for n_prime(0.5), we calculate -x * n(x)
        let expected_n_prime = -x.to_f64().unwrap() * n(x).unwrap().to_f64().unwrap();

        // Compute n_prime(x)
        let calculated_n_prime = n_prime(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n_prime, expected_n_prime, epsilon = 1e-8);
    }

    #[test]
    fn test_n_prime_negative_small_value() {
        // Case where x is a small negative value
        let x = dec!(-0.5);

        // Expected result for n_prime(-0.5), we calculate -x * n(x)
        let expected_n_prime = (-x * n(x).unwrap()).to_f64().unwrap();

        // Compute n_prime(x)
        let calculated_n_prime = n_prime(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n_prime, expected_n_prime, epsilon = 1e-8);
    }

    #[test]
    fn test_n_prime_large_positive_value() {
        // Case where x is a large positive value
        let x = dec!(5.0);

        // Expected result for n_prime(5.0), we calculate -x * n(x)
        let expected_n_prime = -x.to_f64().unwrap() * n(x).unwrap().to_f64().unwrap();

        // Compute n_prime(x)
        let calculated_n_prime = n_prime(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n_prime, expected_n_prime, epsilon = 1e-8);
    }

    #[test]
    fn test_n_prime_large_negative_value() {
        // Case where x is a large negative value
        let x = -dec!(5.0);

        // Expected result for n_prime(-5.0), we calculate -x * n(x)
        let expected_n_prime = (-x * n(x).unwrap()).to_f64().unwrap();

        // Compute n_prime(x)
        let calculated_n_prime = n_prime(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_n_prime, expected_n_prime, epsilon = 1e-8);
    }

    #[test]
    fn test_n_prime_extreme_positive_value() {
        // Case where x is a very large positive value
        let x = dec!(100.0);

        // Expected result for n_prime(100.0), should be extremely close to 0
        let expected_n_prime = (-x * n(x).unwrap()).to_f64().unwrap();

        // Compute n_prime(x)
        let calculated_n_prime = n_prime(x).unwrap().to_f64().unwrap();

        // Assert that n_prime(x) is effectively 0 for such a large input
        assert_relative_eq!(calculated_n_prime, expected_n_prime, epsilon = 1e-100);
    }

    #[test]
    fn test_n_prime_extreme_negative_value() {
        // Case where x is a very large negative value
        let x = -dec!(100.0);

        // Expected result for n_prime(-100.0), should be extremely close to 0
        let expected_n_prime = (-x * n(x).unwrap()).to_f64().unwrap();

        // Compute n_prime(x)
        let calculated_n_prime = n_prime(x).unwrap().to_f64().unwrap();

        // Assert that n_prime(x) is effectively 0 for such a large negative input
        assert_relative_eq!(calculated_n_prime, expected_n_prime, epsilon = 1e-100);
    }
}

#[cfg(test)]
mod calculate_big_n_values {
    use super::*;
    use approx::assert_relative_eq;
    use rust_decimal_macros::dec;
    use statrs::distribution::Normal;

    #[test]
    fn test_big_n_zero() {
        // Case where x = 0.0
        let x = Decimal::ZERO;

        // The CDF of the standard normal distribution at x = 0 is 0.5
        let expected_big_n = 0.5;

        // Compute big_n(x)
        let calculated_big_n = big_n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_big_n, expected_big_n, epsilon = 1e-8);
    }

    #[test]
    fn test_big_n_one() {
        // Case where x = 0.0
        let x = Decimal::ONE;

        // The CDF of the standard normal distribution at x = 1 is 0.841344746054943
        let expected_big_n = 0.841344746054943;

        // Compute big_n(x)
        let calculated_big_n = big_n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_big_n, expected_big_n, epsilon = 1e-8);
    }

    #[test]
    fn test_big_n_two() {
        // Case where x = 0.0
        let x = Decimal::TWO;

        // The CDF of the standard normal distribution at x = 2 is 0.977249868052837
        let expected_big_n = 0.977249868052837;

        // Compute big_n(x)
        let calculated_big_n = big_n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_big_n, expected_big_n, epsilon = 1e-8);
    }

    #[test]
    fn test_big_n_positive_small_value() {
        // Case where x is a small positive value
        let x = dec!(0.5);

        // The expected CDF for the standard normal distribution at x = 0.5 can be precomputed
        let normal_distribution = Normal::new(0.0, 1.0).unwrap();
        let expected_big_n = normal_distribution.cdf(x.to_f64().unwrap());

        // Compute big_n(x)
        let calculated_big_n = big_n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_big_n, expected_big_n, epsilon = 1e-8);
    }

    #[test]
    fn test_big_n_negative_small_value() {
        // Case where x is a small negative value
        let x = -dec!(0.5);

        // The expected CDF for the standard normal distribution at x = -0.5 can be precomputed
        let normal_distribution = Normal::new(0.0, 1.0).unwrap();
        let expected_big_n = normal_distribution.cdf(x.to_f64().unwrap());

        // Compute big_n(x)
        let calculated_big_n = big_n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_big_n, expected_big_n, epsilon = 1e-8);
    }

    #[test]
    fn test_big_n_large_positive_value() {
        // Case where x is a large positive value
        let x = dec!(5.0);

        // The CDF for large positive x should be very close to 1
        let expected_big_n = 1.0f64;

        // Compute big_n(x)
        let calculated_big_n = big_n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_big_n, expected_big_n, epsilon = 1e-6); // if lower epsilon fail
    }

    #[test]
    fn test_big_n_large_negative_value() {
        // Case where x is a large negative value
        let x = -dec!(5.0);

        // The CDF for large negative x should be very close to 0
        let expected_big_n = 0.0f64;

        // Compute big_n(x)
        let calculated_big_n = big_n(x).unwrap().to_f64().unwrap();

        // Assert that the calculated value is close to the expected value
        assert_relative_eq!(calculated_big_n, expected_big_n, epsilon = 1e-6); // if lower epsilon fail
    }

    #[test]
    fn test_big_n_extreme_positive_value() {
        // Case where x is an extremely large positive value
        let x = dec!(100.0);

        // The CDF for an extremely large positive x should be effectively 1
        let expected_big_n = 1.0f64;

        // Compute big_n(x)
        let calculated_big_n = big_n(x).unwrap().to_f64().unwrap();

        // Assert that big_n(x) is effectively 1 for such a large positive input
        assert_relative_eq!(calculated_big_n, expected_big_n, epsilon = 1e-12);
    }

    #[test]
    fn test_big_n_extreme_negative_value() {
        // Case where x is an extremely large negative value
        let x = -dec!(100.0);

        // The CDF for an extremely large negative x should be effectively 0
        let expected_big_n = 0.0f64;

        // Compute big_n(x)
        let calculated_big_n = big_n(x).unwrap().to_f64().unwrap();

        // Assert that big_n(x) is effectively 0 for such a large negative input
        assert_relative_eq!(calculated_big_n, expected_big_n, epsilon = 1e-12);
    }
}

#[cfg(test)]
mod tests_d1_d2_edge_cases {
    use super::*;
    use crate::assert_decimal_eq;
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    #[test]
    fn test_d1_zero_underlying_price() {
        let result = d1(
            Positive::ZERO,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(0.2),
        );
        assert!(result.is_err());
    }

    #[test]
    fn test_d2_negative_rates_and_high_volatility() {
        let result = d2(
            Positive::HUNDRED,
            Positive::HUNDRED,
            -dec!(0.05), // tasa negativa
            Positive::ONE,
            pos_or_panic!(0.8), // alta volatilidad
        )
        .unwrap();
        assert_decimal_eq!(result, dec!(-0.4625), dec!(0.000001));
    }

    #[test]
    fn test_d1_d2_combination_extreme_values() {
        let result_d1 = d1(
            pos_or_panic!(1000.0),
            pos_or_panic!(10.0),
            dec!(0.15),
            Positive::TEN,
            pos_or_panic!(0.9),
        )
        .unwrap();
        let result_d2 = d2(
            pos_or_panic!(1000.0),
            pos_or_panic!(10.0),
            dec!(0.15),
            Positive::TEN,
            pos_or_panic!(0.9),
        )
        .unwrap();
        assert_decimal_eq!(result_d1, dec!(3.5681), dec!(0.0001));
        assert_decimal_eq!(result_d2, dec!(0.7221), dec!(0.0001));
    }
}

#[cfg(test)]
mod tests_probability_density {
    use super::*;
    use rust_decimal_macros::dec;

    #[test]
    fn test_n_prime_symmetry() {
        let x = dec!(1.5);
        let n_prime_pos = n_prime(x).unwrap();
        let n_prime_neg = n_prime(-x).unwrap();
        assert_eq!(n_prime_pos, -n_prime_neg);
    }

    #[test]
    fn test_n_integration_limits() {
        let x_very_large = dec!(10.0);
        let result = n(x_very_large).unwrap();
        assert!(result < dec!(0.0001));
    }

    #[test]
    fn test_n_prime_zero_crossing() {
        let result = n_prime(dec!(0.0)).unwrap();
        assert_eq!(result, dec!(0.0));
    }
}

#[cfg(test)]
mod tests_cumulative_distribution {
    use super::*;
    use rust_decimal_macros::dec;

    #[test]
    fn test_big_n_continuity() {
        let x1 = dec!(0.001);
        let x2 = dec!(-0.001);
        let result1 = big_n(x1).unwrap();
        let result2 = big_n(x2).unwrap();
        assert!((result1 - result2).abs() < dec!(0.001));
    }

    #[test]
    fn test_big_n_conversion_error() {
        let x = Decimal::MAX;
        let result = big_n(x);
        assert!(result.is_ok());
    }

    #[test]
    fn test_big_n_boundary_values() {
        let result_zero = big_n(dec!(0.0)).unwrap();
        assert_eq!(result_zero, dec!(0.5));
    }
}

#[cfg(test)]
mod tests_calculate_d_values_bis {
    use super::*;
    use crate::assert_decimal_eq;
    use crate::model::ExpirationDate;
    use crate::model::types::{OptionStyle, OptionType, Side};
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    #[test]
    fn test_calculate_d_values_with_expiration() {
        let option = Options {
            option_type: OptionType::European,
            side: Side::Long,
            underlying_symbol: "TEST".to_string(),
            strike_price: Positive::HUNDRED,
            underlying_price: Positive::HUNDRED,
            risk_free_rate: dec!(0.05),
            implied_volatility: pos_or_panic!(0.5),
            expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
            quantity: Positive::ONE,
            option_style: OptionStyle::Call,
            dividend_yield: Positive::ZERO,
            exotic_params: None,
        };
        let (d1, d2) = calculate_d_values(&option).unwrap();
        assert_decimal_eq!(d1, dec!(0.1003), dec!(0.0001));
        assert_decimal_eq!(d2, dec!(-0.0430), dec!(0.0001));
    }
}

#[cfg(test)]
mod tests_edge_cases_and_errors {
    use super::*;
    use positive::pos_or_panic;

    use rust_decimal_macros::dec;

    #[test]
    fn test_extreme_volatility_values() {
        let result = d1(
            Positive::HUNDRED,
            Positive::HUNDRED,
            dec!(0.05),
            Positive::ONE,
            pos_or_panic!(1000.0),
        );
        assert!(result.is_ok());
    }

    #[test]
    fn test_precision_limits() {
        let x = dec!(15.0);
        let n_result = n(x).unwrap();
        assert!(n_result.to_f64().unwrap() == 0.0);
    }
}