"""Black-Scholes gamma/delta recompute."""
import numpy as np
from scipy.stats import norm


def bs_gamma(S: float, K: float, T: float, iv: float,
             r: float = 0.04, q: float = 0.0) -> float:
    """Standard Black-Scholes gamma."""
    if T <= 0 or iv <= 0 or S <= 0 or K <= 0:
        return 0.0
    sqrt_T = np.sqrt(T)
    d1 = (np.log(S / K) + (r - q + 0.5 * iv ** 2) * T) / (iv * sqrt_T)
    return np.exp(-q * T) * norm.pdf(d1) / (S * iv * sqrt_T)


def bs_delta(S: float, K: float, T: float, iv: float,
             cp: str, r: float = 0.04, q: float = 0.0) -> float:
    """Black-Scholes delta. cp='C' or 'P'."""
    if T <= 0 or iv <= 0 or S <= 0 or K <= 0:
        return 0.0
    sqrt_T = np.sqrt(T)
    d1 = (np.log(S / K) + (r - q + 0.5 * iv ** 2) * T) / (iv * sqrt_T)
    if cp == "C":
        return np.exp(-q * T) * norm.cdf(d1)
    else:
        return np.exp(-q * T) * (norm.cdf(d1) - 1.0)


def bs_gamma_vec(S_arr, K_arr, T_arr, iv_arr,
                 r: float = 0.04, q: float = 0.0):
    """Vectorised BS gamma over arrays."""
    S_arr = np.asarray(S_arr, dtype=float)
    K_arr = np.asarray(K_arr, dtype=float)
    T_arr = np.asarray(T_arr, dtype=float)
    iv_arr = np.asarray(iv_arr, dtype=float)

    valid = (T_arr > 0) & (iv_arr > 0) & (S_arr > 0) & (K_arr > 0)
    out = np.zeros_like(S_arr)
    if valid.any():
        sqrt_T = np.sqrt(T_arr[valid])
        d1 = (np.log(S_arr[valid] / K_arr[valid]) +
              (r - q + 0.5 * iv_arr[valid] ** 2) * T_arr[valid]) / (iv_arr[valid] * sqrt_T)
        out[valid] = np.exp(-q * T_arr[valid]) * norm.pdf(d1) / (S_arr[valid] * iv_arr[valid] * sqrt_T)
    return out


def bs_delta_2d(s, K, T, iv, cp_sign, r: float = 0.04, q: float = 0.0):
    """Vectorised BS delta over a 2D broadcast.

    s      : (1, G) candidate spot levels
    K,T,iv : (n, 1) per-contract arrays
    cp_sign: (n, 1) +1 for calls, -1 for puts
    Returns (n, G) deltas. Call delta positive, put delta negative (no extra sign).
    """
    sqrt_T = np.sqrt(T)
    d1 = (np.log(s / K) + (r - q + 0.5 * iv ** 2) * T) / (iv * sqrt_T)
    # call delta = e^{-qT} N(d1); put delta = e^{-qT} (N(d1) - 1)
    # combine via cp_sign: delta = e^{-qT} * ( N(d1) - (cp_sign<0) )
    is_put = (cp_sign < 0)
    return np.exp(-q * T) * (norm.cdf(d1) - is_put)