"""European options under Black-Scholes, with a continuous dividend yield.

Mutants are made by editing this source text, so each line that a mutant
replaces has to stay unique. See mutants.py.
"""
from math import exp, log, pi, sqrt
from statistics import NormalDist

_normal = NormalDist()
N = _normal.cdf


def n(x):
    return exp(-0.5 * x * x) / sqrt(2.0 * pi)


def d1_d2(spot, strike, rate, div, vol, time):
    sigma_root = vol * sqrt(time)
    d1 = (log(spot / strike) + (rate - div + 0.5 * vol * vol) * time) / sigma_root
    d2 = d1 - vol * sqrt(time)
    return d1, d2


def forward_value(spot, strike, rate, div, time):
    return spot * exp(-div * time), strike * exp(-rate * time)


def is_degenerate(vol, time):
    return vol * sqrt(time) < 1e-12


def call(spot, strike, rate, div, vol, time):
    carried, discounted = forward_value(spot, strike, rate, div, time)
    if is_degenerate(vol, time):
        return max(carried - discounted, 0.0)
    d1, d2 = d1_d2(spot, strike, rate, div, vol, time)
    return carried * N(d1) - discounted * N(d2)


def put(spot, strike, rate, div, vol, time):
    carried, discounted = forward_value(spot, strike, rate, div, time)
    if is_degenerate(vol, time):
        return max(discounted - carried, 0.0)
    d1, d2 = d1_d2(spot, strike, rate, div, vol, time)
    return discounted * N(-d2) - carried * N(-d1)


def delta(spot, strike, rate, div, vol, time):
    d1, _ = d1_d2(spot, strike, rate, div, vol, time)
    return exp(-div * time) * N(d1)


def gamma(spot, strike, rate, div, vol, time):
    d1, _ = d1_d2(spot, strike, rate, div, vol, time)
    return exp(-div * time) * n(d1) / (spot * vol * sqrt(time))


def vega(spot, strike, rate, div, vol, time):
    d1, _ = d1_d2(spot, strike, rate, div, vol, time)
    return spot * exp(-div * time) * n(d1) * sqrt(time)


def implied_vol(target, spot, strike, rate, div, time):
    """Bisection on the call price. Returns None if the target is out of range."""
    low, high = 1e-6, 5.0
    if not call(spot, strike, rate, div, low, time) <= target <= call(spot, strike, rate, div, high, time):
        return None
    for _ in range(200):
        middle = 0.5 * (low + high)
        if call(spot, strike, rate, div, middle, time) < target:
            low = middle
        else:
            high = middle
    return 0.5 * (low + high)