//
// This source code resides at www.jaeckel.org/StrikeFromVolatilityAndDeltaWithPremium.7z .
//
// ======================================================================================
// Copyright © 2020 Peter Jäckel.
// 
// Permission to use, copy, modify, and distribute this software is freely granted,
// provided that this notice is preserved.
//
// WARRANTY DISCLAIMER
// The Software is provided "as is" without warranty of any kind, either express or implied,
// including without limitation any implied warranties of condition, uninterrupted use,
// merchantability, fitness for a particular purpose, or non-infringement.
// ======================================================================================
//

#define NOMINMAX // to suppress MSVC's definition of min() and max()
#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include <cmath>
#include <limits>
#include <algorithm>
#include <assert.h>

#include "normaldistribution.h"
#include "lets_be_rational.h"

#include "XLFunctions.h"

#ifdef _MSC_VER
# define snprintf(buffer,sizeOfBuffer,...) _snprintf_s((buffer),(sizeOfBuffer),_TRUNCATE,__VA_ARGS__)
#endif

extern "C" double InverseModifiedDeltaWithPremiumFunction(double f, double alpha);

double AttainableAccuracyFactor(double y, double alpha);

#define SIGN( X ) ( X < 0 ? -1 : ( X > 0 ? 1 : 0 ) )

// See http://gcc.gnu.org/wiki/Visibility
#ifndef EXPORT_EXTERN_C
# if defined( WIN32 ) || defined( _WIN32 ) || defined( WIN64 ) || defined( _WIN64 )
#  define EXPORT_EXTERN_C extern "C" __declspec(dllexport)
# else
#  define EXPORT_EXTERN_C extern "C" __attribute__((visibility("default")))
# endif
#endif

#define LN_TWO    0.69314718055994530941723212145817656807550013436026
#define LN_TWO_PI 1.8378770664093453390819377091247588396072387695313
#define SQRT_PI   1.77245385090551602729816748334
#define SQRT_TWO         1.41421356237309504880168872421
#define SQRT_TWO_OVER_PI 0.7978845608028653558798921198687637

namespace {
    const double one_over_sqrt_dbl_epsilon = 1/sqrt(DBL_EPSILON);
    const double NaN = std::numeric_limits<double>::quiet_NaN();
}

DECLARE_XL_FUNCTION( norm_cdf,         ".NormsDist" , "BB", "x",                      "the cumulative normal probability distribution function.", "the argument x to Phi(x).")
DECLARE_XL_FUNCTION( inverse_norm_cdf, ".NormsInv" ,  "BB", "p",                      "the inverse cumulative normal probability distribution function.", "the probability p such that the returned value x matches p=Phi(x).")
DECLARE_XL_FUNCTION( erfc_cody,        ".erfc",       "BB", "x",                      "the complementary error function erfc().", "the argument x to erfc(x).")
DECLARE_XL_FUNCTION( erfcx_cody,       ".erfcx",      "BB", "x",                      "the scaled complementary error function erfcx().", "the argument x to erfcx(x).")
DECLARE_XL_FUNCTION( black,            ".Black",  "BBBBBB", "F,K,volatility,T,theta", "returns the Black (forward) option price.", "the forward", "the strike", "the volatility", "the time to expiry", "+/-1 for calls/puts")

EXPORT_EXTERN_C double SqrtOnePlusXMinusOne(double x) {
    if ( fabs (x) < 0.25) { // sqrt(1+x) -1 expanded for small x. PJ , 2019.
        const double g = x < 0
            ? // net relative accuracy [of x*(0.5-x*g)] better than 9E-17
            (1.249999999999904E-1+x*(1.952674587134858E-1+x*(8.609312134216375E-2+9.372014538795577E-3*x)))/(1+x*(2.06213966971922+x*(1.407314806167365+x*(3.529648863602843E-1+2.372729847070329E-2*x))))
            : // net relative accuracy [of x*(0.5-x*g)] better than 8E-18
        (1.24999999999998E-1+x*(1.811342940763106E-1+x*(7.189411491089435E-2+6.684805698482948E-3*x)))/(1+x*(1.949074352607691+x*(1.237190095758854+x*(2.817377533594094E-1+1.654456755072703E-2*x))));
        return x * ( 0.5 - x*g );
    }
    return sqrt(1.0 + x) - 1.0;
}

DECLARE_XL_FUNCTION( SqrtOnePlusXMinusOne, ".SqrtOnePlusXMinusOne",  "BB", "x", "sqrt(1+x)-1 with accurate rational function approximation when -0.25<x<0.25.", "the argument x in sqrt(1+x)-1.")

// From the source code of "Let's Be Rational".
double AsymptoticExpansionOfCdfOverPdf(double z) {
    // The asymptotic expansion  Φ(z) = φ(z)/|z|·[1-1/z^2+...],  Abramowitz & Stegun (26.2.12), suffices for Φ(z) to have

    // Asymptotic expansion for very negative z following (26.2.12) on page 408
    // in M. Abramowitz and A. Stegun, Pocketbook of Mathematical Functions, ISBN 3-87144818-4, is
    //   Φ(z) = φ(z)/(-z)·[1-1/z^2+...].
    // Hence, 
    //   Φ(z) / φ(z) = -(1-1/z^2+...) / z.
    double sum = 1;
    if (z >= -one_over_sqrt_dbl_epsilon) {
        double zsqr = z * z, i = 1, g = 1, x, y, a = DBL_MAX, lasta;
        do {
            lasta = a;
            x = (4 * i - 3) / zsqr;
            y = x * ((4 * i - 1) / zsqr);
            a = g * (x - y);
            sum -= a;
            g *= y;
            ++i;
            a = fabs(a);
        } while (lasta > a && a >= fabs(sum * DBL_EPSILON));
    }
    return -sum / z;
}

EXPORT_EXTERN_C double LnOfOnePlusX(double x) {
    if (fabs(x) < 0.25){
        // Nonlinear-Remez optimised Chebyshev-Padé(5,5) of theoretical relative accuracy better than 5E-18.
        const double g = 1 - x * ((4.9999999999999978537E-1 + x * (1.1365571796459058354E0 + x * (8.8850223817666576914E-1 + x * (2.741489129503292182E-1 + x * (2.6720133887436275504E-2 + 3.3286545792509327678E-5 * x))))) / (1 + x * (2.9397810259584798844E0 + x * (3.2368584936588268311E0 + x * (1.6363129753602314485E0 + x * (3.6846541491884293455E-1 + 2.8084369911925627839E-2 * x))))));
        return x * g;
    }
    return log(1 + x);
}

EXPORT_EXTERN_C double LnCdf(double z) { // ln(Φ(z))
    if (z > 0)
        return LnOfOnePlusX(-norm_cdf(-z));
    if (z <= -9.0)
        return log(AsymptoticExpansionOfCdfOverPdf(z)) - 0.5 * (z * z + LN_TWO_PI);
    return log(norm_cdf(z));
}

EXPORT_EXTERN_C double LnTwoCdf(double y) { // ln(2·Φ(y))
    // For |y|<0.4, we use a Nonlinear-Remez optimized approximation given by f(y)=θ(y,g(y)) with θ(y,g) = y·(√(2/π)-y·g) and the (5,4)-rational function
    //   g(y) = (3.183098861837894E-1+y*(-2.1883026451258495E-1+y*(7.2342988529039326E-2+y*(-1.1919277649580325E-2+y*(8.3837904820976024E-4-1.2552162258772677E-6*y)))))/(1+y*(-5.7332388994577462E-1+y*(1.768497738757971E-1+y*(-2.598729916076733E-2+1.7247546759678694E-3*y))))
    // The net accuracy is better than 1.122E-16.
    return fabs(y) > 0.4 ? LN_TWO + LnCdf(y) : y * (0.79788456080286535588 - y * ((3.183098861837894E-1 + y * (-2.1883026451258495E-1 + y * (7.2342988529039326E-2 + y * (-1.1919277649580325E-2 + y * (8.3837904820976024E-4 - 1.2552162258772677E-6 * y))))) / (1 + y * (-5.7332388994577462E-1 + y * (1.768497738757971E-1 + y * (-2.598729916076733E-2 + 1.7247546759678694E-3 * y))))));
}

EXPORT_EXTERN_C double InverseLnErfcx(double y) { // Inverse of y(z) := ln(erfcx(z))
    double z;
    //
    // All of the below rational functions are Remez-II(b) optimized on their respective intervals.
    //
    if (y >= 2.926054502) {
        // z ∈ [-∞,-3/2] <══> y ∈ [2.926054502,∞] <══> ω ∈ [0,0.4478466214].
        // z(y) = θ(ω,ζ) := -(1+exp(-ζ/ω)/(4√π))/√ω with ω:=1/(y-ln(2)) and ζ=ζ(ω) being a (7,7)-rational function.
        // Relative accuracy of z(y): 2.1E-17 (in perfect arithmetic).
        // Since GCC [7.4] does not permit unicode identifiers such as ω, we use w instead of ω, and zeta instead of ζ.
        const double w = 1 / (y - LN_TWO);
        const double zeta = (1.008831813880055711E0 + w * (5.6930462040180993295E1 + w * (7.8758967157473104837E2 + w * (6.5369557790922766761E3 + w * (3.207323236990434636E4 + w * (4.7081311888818542321E4 + w * ((-1.5356447406703549185E3) - 1.3618612709837885034E4 * w))))))) / (1 + w * (4.9386339997921838307E1 + w * (5.6659651701221432338E2 + w * (4.6794735596030363518E3 + w * (1.8438296226659342414E4 + w * (2.3066435716227694152E4 + w * (6.0425284640091166861E3 - 1.8361146033364165539E2 * w)))))));
        z = -(1 + exp(-zeta / w) / (4 * SQRT_PI)) / sqrt(w);
    } else if (y > -1.134492089) {
        // z ∈ [-3/2,3/2] <══> y ∈ [-1.134492089,2.926054502]
        // z(y) = θ(y,γ) := y·(z_1+y·γ) with γ=γ(y) being a (9,9)-rational function and z_1 := -(√π)/2.
        // Relative accuracy of z(y): 6.2E-18 (in perfect arithmetic).
        // Since GCC [7.4] does not permit unicode identifiers, we use gamma instead of γ.
        const double gamma = (2.5292753687758446586E-1 + y * (3.215776976543350707E-1 + y * (2.1543770801630490759E-1 + y * (8.8809464227384136275E-2 + y * (2.473545115937593774E-2 + y * (4.6565738088983181102E-3 + y * (5.9311148077307959971E-4 + y * (4.6804490045539737994E-5 + y * (1.952362913447388075E-6 + 1.3547814811265750855E-9 * y))))))))) / (1 + y * (1.5915724440392176382E0 + y * (1.2703842021304194846E0 + y * (6.3481645737589602599E-1 + y * (2.1605555010263235345E-1 + y * (5.1491662689748139232E-2 + y * (8.5935378171905465074E-3 + y * (9.7598654688020152114E-4 + y * (6.9714835279954410095E-5 + 2.4561708908515460842E-6 * y)))))))));
        z = y * (-0.5 * SQRT_PI + y * gamma);
    } else {
        // z ∈ [3/2,+∞] <══> y ∈ [-∞,-1.134492089]
        // z(y) = θ(μ,R) := (1-μ·R)/√(μ·π) with R=R(μ) being a (7,7)-rational function, and μ := exp(2·y).
        // Relative accuracy of z(y) is 1.01E-17 (in perfect arithmetic).
        // Since GCC [7.4] does not permit unicode identifiers, we use m instead of μ.
        const double m = exp(2 * y);
        const double R = (1.5707963267948264287E0 + m * (1.4881863687985725824E2 + m * (5.0573969043451462837E3 + m * (7.71433208578163961E4 + m * (5.4104053714173542146E5 + m * (1.5996711891638245768E6 + m * (1.5582800959832326392E6 + 2.2979729179982853295E5 * m))))))) / (1 + m * (9.6311683061208178375E1 + m * (3.3610552996661054666E3 + m * (5.3544584629105503552E4 + m * (4.039412696733845392E5 + m * (1.3577968131365781872E6 + m * (1.6990431126157576341E6 + 4.9222272232082468876E5 * m)))))));
        z = (1 - m * R) / (SQRT_PI * sqrt(m));
    }
    return z;
}

EXPORT_EXTERN_C double ForwardDelta(double F, double K, double sigma, double T, double q /* q=±1 */, short /* Excel's BOOL type (code 'A') is a "short" */ with_premium) {
    const double s = fabs(sigma) * sqrt(T);
    if (s > DBL_MIN) {
        const double delta = q * norm_cdf((F > 0 ? q : -q) * (log(F / K) / s + 0.5 * s));
        return with_premium ? delta - black(F,K,sigma,T,q)/F : delta;
    }
    const double intrinsic = std::max(SIGN(q)*(F-K),0.0);
    const double delta = ( intrinsic > 0 ? SIGN(q) : 0 );
    return with_premium ? delta - intrinsic/F : delta;
}

DECLARE_XL_FUNCTION( ForwardDelta, ".ForwardDelta",  "BBBBBBA", "F,K,volatility,T,theta,with_premium", "returns the Black (forward) delta.", "the forward", "the strike", "the volatility", "the time to expiry", "+/-1 for calls/puts", "TRUE if the delta is to be with premium, else FALSE")

EXPORT_EXTERN_C double MinimumCallStrikeForForwardDeltaWithPremium(double volatility, double F, double T) {
    const double stddev = volatility * sqrt(T);
    const double ymax = SQRT_TWO * InverseLnErfcx(log(SQRT_TWO_OVER_PI / stddev));
    return F * exp(stddev * (ymax - 0.5 * stddev));
}

DECLARE_XL_FUNCTION( MinimumCallStrikeForForwardDeltaWithPremium, ".MinimumCallStrikeForForwardDeltaWithPremium",  "BBBB", "volatility,F,T", "returns the lowest call option strike for which a forward delta with premium can be resolved to the same strike.", "the volatility", "the forward", "the time to expiry")

EXPORT_EXTERN_C double AttainableAccuracy(double F, double K, double volatility, double T, double q /* q=±1 */) {
    const double alpha = SIGN( q ) * volatility * sqrt(T),y = log(K/F) / alpha + 0.5 * alpha;
    return AttainableAccuracyFactor(y,alpha) * DBL_EPSILON;
}

DECLARE_XL_FUNCTION( AttainableAccuracy, ".AttainableAccuracy",  "BBBBBB", "F,K,volatility,T,theta", "returns an estimate for the numerically attainable relative accuracy given by max( |f(y)/(y·f'(y)| , 1 ).", "the forward", "the strike", "the volatility", "the time to expiry", "+/-1 for calls/puts")

double StrikeFromForwardDeltaWithoutPremium(double volatility, double forward_delta, double F, double T) {
    const double ad = fabs(forward_delta);
    if (!(ad > 0 && ad < 1))
        return NaN;
    return F * exp(volatility * (0.5 * volatility * T - SIGN( forward_delta ) * sqrt(T) * inverse_norm_cdf(ad)));
}

EXPORT_EXTERN_C double StrikeFromForwardDelta(double volatility, double forward_delta, double F, double T, short /* Excel's BOOL type (code 'A') is a "short" */ with_premium) {
    if (with_premium) {
        // See "Strike from volatility and delta with premium" for details.
        if (!(forward_delta < 1))
            return NaN;
        const double alpha = SIGN( forward_delta ) * volatility * sqrt(T), half_alpha_square = 0.5 * alpha * alpha;
        const double y = InverseModifiedDeltaWithPremiumFunction(log(2 * fabs(forward_delta)) + half_alpha_square, alpha);
        const double strike = F * exp(alpha * y - half_alpha_square);
        return strike;
    }
    return StrikeFromForwardDeltaWithoutPremium(volatility, forward_delta, F, T);
}

const char *StrikeFromForwardDeltaSynonyms[] = { ".StrikeFromForwardDelta", "StrikeFromForwardDelta", 0 };

DECLARE_XL_FUNCTION_WITH_SYNOYMS( StrikeFromForwardDelta, StrikeFromForwardDeltaSynonyms, "BBBBBA", "volatility,forward_delta,F,T,with_premium", "returns the strike K from a forward delta, volatility, the forward F, time to expiry T, where the delta is either 'with premium' or without.", "the volatility.", "the forward delta", "the forward", "the time to expiry", "TRUE if the delta is with premium, else FALSE")
