import numpy as np
from scipy.stats import norm


def black_scholes(
    S: float,
    K: float,
    T: float,
    r: float,
    sigma: float,
    option_type: str = "call",
) -> dict:
    """
    Black-Scholes European option pricer.

    Parameters
    ----------
    S           : float — Spot price          (e.g. 100.0)
    K           : float — Strike price         (e.g. 105.0)
    T           : float — Time to expiry, years (e.g. 0.5)
    r           : float — Risk-free rate p.a.  (e.g. 0.05)
    sigma       : float — Implied vol p.a.     (e.g. 0.20)
    option_type : str   — "call" or "put"

    Returns
    -------
    dict with: price, delta, gamma, vega, theta, rho, d1, d2
    """
    if T <= 0:
        intrinsic = (
            max(S - K, 0.0) if option_type == "call" else max(K - S, 0.0)
        )
        return {
            "price": intrinsic,
            "delta": 1.0 if (option_type == "call" and S > K) else (
                     -1.0 if (option_type == "put" and S < K) else 0.0),
            "gamma": 0.0, "vega": 0.0, "theta": 0.0, "rho": 0.0,
            "d1": None, "d2": None,
        }

    sqrtT = np.sqrt(T)
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sqrtT)
    d2 = d1 - sigma * sqrtT

    Nd1  = norm.cdf(d1)   # N(d1)  — risk-neutral delta of call
    Nd2  = norm.cdf(d2)   # N(d2)  — RN probability of finishing ITM
    nd1  = norm.pdf(d1)   # φ(d1)  — standard normal PDF
    disc = np.exp(-r * T) # e^{-rT} — discount factor

    if option_type == "call":
        price = S * Nd1 - K * disc * Nd2
        delta = Nd1
        rho   = K * T * disc * Nd2 / 100        # per 1 bp move in r
    else:
        price = K * disc * (1 - Nd2) - S * (1 - Nd1)
        delta = Nd1 - 1
        rho   = -K * T * disc * (1 - Nd2) / 100

    gamma = nd1 / (S * sigma * sqrtT)
    vega  = S * nd1 * sqrtT / 100               # per 1% move in sigma
    theta = (
        -S * nd1 * sigma / (2 * sqrtT)
        - r * K * disc * (Nd2 if option_type == "call" else (1 - Nd2))
    ) / 365                                      # per calendar day

    return {
        "price": round(price, 4),
        "delta": round(delta, 4),
        "gamma": round(gamma, 6),
        "vega":  round(vega,  4),
        "theta": round(theta, 4),
        "rho":   round(rho,   4),
        "d1":    round(d1, 4),
        "d2":    round(d2, 4),
    }


def implied_vol(
    market_price: float,
    S: float,
    K: float,
    T: float,
    r: float,
    option_type: str = "call",
    tol: float = 1e-6,
    max_iter: int = 200,
) -> float | None:
    """Newton-Raphson inversion to find implied volatility."""
    sigma = 0.20  # initial guess
    for _ in range(max_iter):
        result = black_scholes(S, K, T, r, sigma, option_type)
        diff   = result["price"] - market_price
        vega   = result["vega"] * 100          # convert back to raw vega
        if abs(vega) < 1e-10:
            break
        sigma -= diff / vega
        sigma  = max(sigma, 1e-6)              # keep sigma positive
        if abs(diff) < tol:
            return round(sigma, 6)
    return None


# ── Example usage ───────────────────────────────────────────────
if __name__ == "__main__":
    # Scenario 1 : ATM call — 6-month, AAPL-style
    atm = black_scholes(S=180, K=180, T=0.5, r=0.05, sigma=0.25)
    print("=== ATM Call (S=K=180, T=6m, r=5%, σ=25%) ===")
    print(f"  Price : {atm['price']:.4f}")
    print(f"  Delta : {atm['delta']:.4f}  (≈ 0.50 pour ATM)")
    print(f"  Gamma : {atm['gamma']:.6f}")
    print(f"  Vega  : {atm['vega']:.4f}  (per 1% σ)")
    print(f"  Theta : {atm['theta']:.4f}  (per day)")
    print()

    # Scenario 2 : Put-Call Parity check
    call = black_scholes(S=180, K=180, T=0.5, r=0.05, sigma=0.25)["price"]
    put  = black_scholes(S=180, K=180, T=0.5, r=0.05, sigma=0.25,
                         option_type="put")["price"]
    parity = 180 - 180 * np.exp(-0.05 * 0.5)
    print("=== Put-Call Parity ===")
    print(f"  Call - Put      : {call - put:.4f}")
    print(f"  S - K·e^(-rT)  : {parity:.4f}")
    print(f"  Parity holds   : {abs((call - put) - parity) < 1e-8}")
    print()

    # Scenario 3 : Implied vol inversion
    iv = implied_vol(market_price=18.50, S=180, K=180, T=0.5, r=0.05)
    print(f"=== Implied Vol ===")
    print(f"  Market price 18.50 → σ_impl = {iv:.4f} ({iv*100:.2f}%)")

#include <cmath>
#include <string>
#include <stdexcept>
#include <iostream>
#include <iomanip>

// ── Normal distribution helpers ─────────────────────────────────
// Hart (1968) rational approximation — max error 7.5e-8
double norm_cdf(double x) {
    static const double a[] = {
        0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429
    };
    double k    = 1.0 / (1.0 + 0.2316419 * std::abs(x));
    double poly = k * (a[0] + k * (a[1] + k * (a[2] + k * (a[3] + k * a[4]))));
    double phi  = std::exp(-0.5 * x * x) / std::sqrt(2.0 * M_PI);
    double res  = 1.0 - phi * poly;
    return x >= 0.0 ? res : 1.0 - res;
}

double norm_pdf(double x) {
    return std::exp(-0.5 * x * x) / std::sqrt(2.0 * M_PI);
}

// ── Greeks struct ───────────────────────────────────────────────
struct BSResult {
    double price, delta, gamma, vega, theta, rho, d1, d2;
};

// ── Black-Scholes pricer ────────────────────────────────────────
BSResult black_scholes(
    double S, double K, double T, double r, double sigma,
    const std::string& type = "call"
) {
    if (T <= 0.0) {
        double intrinsic = (type == "call")
            ? std::max(S - K, 0.0) : std::max(K - S, 0.0);
        return {intrinsic, (S > K ? 1.0 : 0.0), 0, 0, 0, 0, 0, 0};
    }

    double sqrtT = std::sqrt(T);
    double d1    = (std::log(S / K) + (r + 0.5 * sigma * sigma) * T)
                   / (sigma * sqrtT);
    double d2    = d1 - sigma * sqrtT;
    double disc  = std::exp(-r * T);
    double nd1   = norm_pdf(d1);
    double Nd1   = norm_cdf(d1);
    double Nd2   = norm_cdf(d2);
    bool   isC   = (type == "call");

    double price  = isC ? (S * Nd1 - K * disc * Nd2)
                        : (K * disc * (1.0 - Nd2) - S * (1.0 - Nd1));
    double delta  = isC ? Nd1 : Nd1 - 1.0;
    double gamma  = nd1 / (S * sigma * sqrtT);
    double vega   = S * nd1 * sqrtT / 100.0;
    double theta  = (-S * nd1 * sigma / (2.0 * sqrtT)
                    - r * K * disc * (isC ? Nd2 : (1.0 - Nd2))) / 365.0;
    double rho    = (isC ? K * T * disc * Nd2
                         : -K * T * disc * (1.0 - Nd2)) / 100.0;

    return {price, delta, gamma, vega, theta, rho, d1, d2};
}

// ── Newton-Raphson implied vol ──────────────────────────────────
double implied_vol(
    double market_price, double S, double K, double T, double r,
    const std::string& type = "call",
    double tol = 1e-6, int max_iter = 200
) {
    double sigma = 0.20;
    for (int i = 0; i < max_iter; ++i) {
        auto   res  = black_scholes(S, K, T, r, sigma, type);
        double diff = res.price - market_price;
        double vega = res.vega * 100.0;
        if (std::abs(vega) < 1e-10) break;
        sigma -= diff / vega;
        sigma  = std::max(sigma, 1e-6);
        if (std::abs(diff) < tol) return sigma;
    }
    return -1.0; // failed to converge
}

// ── Main ────────────────────────────────────────────────────────
int main() {
    auto res = black_scholes(180.0, 180.0, 0.5, 0.05, 0.25, "call");
    std::cout << std::fixed << std::setprecision(4);
    std::cout << "=== ATM Call (S=K=180, T=6m, r=5%, sigma=25%) ===\n";
    std::cout << "  Price : " << res.price << "\n";
    std::cout << "  Delta : " << res.delta << "\n";
    std::cout << "  Gamma : " << std::setprecision(6) << res.gamma << "\n";
    std::cout << "  Vega  : " << std::setprecision(4) << res.vega  << "\n";
    std::cout << "  Theta : " << res.theta << "\n\n";

    double iv = implied_vol(18.50, 180.0, 180.0, 0.5, 0.05);
    std::cout << "=== Implied Vol ===\n";
    std::cout << "  Price 18.50 -> sigma_impl = "
              << std::setprecision(6) << iv << "\n";
    return 0;
}