//
// 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 <float.h>
#include <math.h>
#include <cmath>
#include <algorithm>
#include <limits>
#include <assert.h>

#include "normaldistribution.h"
#include "rationalcubic.h"

#include "XLFunctions.h"

// 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

extern "C" double LnTwoCdf(double y);
extern "C" double SqrtOnePlusXMinusOne(double x);
extern "C" double InverseLnErfcx(double y);

#define LN_TWO           0.69314718055994530941723212145817656807550013436026
#define SQRT_TWO         1.41421356237309504880168872421
#define SQRT_TWO_OVER_PI 0.7978845608028653558798921198687637

namespace {
    const double NaN = std::numeric_limits<double>::quiet_NaN();
}

EXPORT_EXTERN_C double HouseHolder6Step(double nu, double h2, double h3, double h4, double h5, double h6) { return (nu * (1 + nu * (2 * h2 + nu * ((3 * (h2 * h2)) / 4 + h3 / 2 + nu * ((h2 * h3) / 6 + h4 / 12 + (h5 * nu) / 120))))) / (1 + nu * ((5 * h2) / 2 + nu * ((3 * (h2 * h2)) / 2 + (2 * h3) / 3 + nu * ((h2 * h2 * h2) / 8 + (h2 * h3) / 2 + h4 / 8 + nu * ((h3 * h3) / 36 + (h2 * h4) / 24 + h5 / 60 + (h6 * nu) / 720))))); }

double is_too_small(double x) { return fabs(x) < DBL_MIN; }

namespace InterpolationRationalCubic {
    const double r_min = -(1 - sqrt(DBL_EPSILON));
    const double r_max = 2 / (DBL_EPSILON * DBL_EPSILON);
    double ControlParameterToFitInteriorPoint(double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double x_interior, double y_interior) {
        assert(x_interior >= x_l && x_r >= x_interior);
        const double h = (x_r - x_l);
        if (is_too_small(h))
            return 3;
        const double t = (x_interior - x_l) / h, omt = 1 - t, t2 = t * t, omt2 = omt * omt, y_abs_max = std::max(fabs(y_interior), std::max(fabs(y_l), fabs(y_r)));
        const double cubic_form = y_r * t2 * t + (3 * y_r - h * d_r) * t2 * omt + (3 * y_l + h * d_l) * t * omt2 + y_l * omt2 * omt, y_interior_minus_cubic_form = y_interior - cubic_form;
        if (fabs(y_interior_minus_cubic_form) <= DBL_EPSILON * y_abs_max)
            return 3;
        const double linear_form = y_r * t + y_l * omt, y_interior_minus_linear_form = y_interior - linear_form;
        if (fabs(y_interior_minus_linear_form) <= DBL_EPSILON * y_abs_max)
            return r_max;
        const double w = t * omt;
        if (is_too_small(w))
            return 3;
        // Maxima code to confirm the below simplified formula for the solution for r:
        //     cubic(r) := y_r * t*t * t + (r * y_r - h * d_r) * t*t * omt + (r * y_l + h * d_l) * t * omt*omt + y_l * omt*omt * omt; linear : y_r*t+y_l*omt;
        //     r : ev(r,solve(cubic(r) - y_interior*(1 + (r - 3) * t * omt),r)[1]); simplified_r : fullratsimp(3-((y_interior-cubic(3))/(y_interior-linear))/(t*omt)); fullratsimp(r-simplified_r);
        double r = 3 - (y_interior_minus_cubic_form / y_interior_minus_linear_form) / w;
        r = std::min(r_max, std::max(r_min, r));
        return r;
    }
}


namespace /* ModifiedDeltaWithPremium */ { // f(y) = ln(2·Φ(-y)) + α·y  as defined in equation (1.14)
    const unsigned int NY = 11;
    const double Yeps = -8.2535; // ≈ Φ^(-1)(DBL_EPSILON·ln(2)/2)
    const double Y[NY] = { -8.2535, -7.0, -6.3333333333333333333, -5.6666666666666666667, -5, -4.3333333333333333333, -3.6666666666666666667, -3.0, -2.0, -1.0, 0.0 };
    // c  :  f(y[i],alpha)=c[i]+y[i]*α
    const double C[NY] = { 6.93147180559945233E-1, 6.9314718055866549687E-1, 6.9314718043998508325E-1, 6.9314717327983520900E-1, 6.9314689390833234565E-1, 6.931398371091453459E-1, 6.9302430662128693903E-1, 6.9179637059519711562E-1, 6.7013427123098182E-1, 5.20393401536495419E-1, 0.0 };
    // d  :  f'(y[i];α) = d[i]+α
    const double D[NY] = { -6.4384757653278754E-16, -9.1347204083762840726E-12, -7.7784480690027837163E-10, -4.2468887773727108186E-8, -1.4867199409049057124E-6, -3.3371107453824317093E-5, -4.8032966799311978792E-4, -4.4378390421256637933E-3, -5.52478626789899591E-2, -2.87599970939178361E-1, -7.97884560802865355E-1 };
    // e: f''(y[i];α) = e[i] (not dependent on α)
    const double E[NY] = { -5.31399597291336238E-15, -6.3943042858717431625E-11, -4.926350444306805564E-9, -2.4065703252139337512E-7, -7.4336019148607112465E-6, -1.4460924593071806876E-4, -1.7614394992313936031E-3, -1.3333211541740806209E-2, -1.13548051688576449E-1, -3.70313714223394599E-1, -6.36619772367581343E-1 };
    // α such that ymax=y : alpha_of_ymax(y)=sqrt(2.0/pi)/erfcx(y/sqrt(2.0))   ---   UNUSED !
    const double αmax[NY] = { 6.4384757653278754E-16, 9.1347204083762840726E-12, 7.7784480690027837163E-10, 4.2468887773727108186E-8, 1.4867199409049057124E-6, 3.3371107453824317093E-5, 4.8032966799311978792E-4, 4.4378390421256637933E-3, 5.5247862678989959102E-2, 2.8759997093917836123E-1, 7.9788456080286535588E-1 };
    double F(unsigned int i, double α) { return C[i] + Y[i] * α; }
    double Fp(unsigned int i, double α) { return D[i] + α; }
    // g[i] = y[i]-1/3
    double G[NY] = { -8.58683333333333333, -7.33333333333333333, -6.66666666666666667, -6.0, -5.33333333333333333, -4.66666666666666667, -4.0, -3.33333333333333333, -2.33333333333333333, -1.33333333333333333, -3.33333333333333333E-1 };
    // h  :  f(g[i],alpha)=h[i]+g[i]*α
    double H[NY] = { 6.93147180559945305E-1, 6.93147180559833061E-1, 6.93147180546861385E-1, 6.93147179573357664E-1, 6.93147132346910496E-1, 6.93145649932037368E-1, 6.9311550881656782E-1, 6.92718028154026267E-1, 6.8328336404937177E-1, 5.97504603819424651E-1, 2.31998089638832187E-1 };
    // Fi(i,α) = F(y[i]-1/3,α) = F(g[i],α)
    double Fi(unsigned int i, double α) { return H[i] + G[i] * α; }
    unsigned int FindInterval(double f, double α, unsigned int i_l, unsigned int i_r, double ymax, double fmax, double& f_y_l, double& f_y_r) {
        // 'left' and 'right' here refers to the ordering in y.
        const unsigned int i_mid = (i_l + i_r) / 2;
        // Recall that f(y) is (supposed to be) downward sloping.
        const double f_y_mid = Y[i_mid] < ymax ? fmax : C[i_mid] + Y[i_mid] * α;
        if (f > f_y_mid) { // branch to the left in y to get the upper branch in f
            f_y_r = f_y_mid;
            return i_l + 1 >= i_mid ? i_l : FindInterval(f, α, i_l, i_mid, ymax, fmax, f_y_l, f_y_r);
        }
        // branch to the right in y to get the lower branch in f
        f_y_l = f_y_mid;
        return i_mid + 1 >= i_r ? i_mid : FindInterval(f, α, i_mid, i_r, ymax, fmax, f_y_l, f_y_r);
    }
    unsigned int FindInterval(double f, double α, double ymax, double fmax, double& f_y_l, double& f_y_r) {
        return FindInterval(f, α, 0, NY, ymax, fmax, f_y_l, f_y_r);
    }
    unsigned int FindYIndexNotLess(double y, unsigned int i_l, unsigned int i_r) {
        if (i_l >= i_r) return i_l;
        const unsigned int i_mid = (i_l + i_r) / 2;
        return y > Y[i_mid] ? (i_mid + 1 >= i_r ? i_r : FindYIndexNotLess(y, i_mid, i_r)) : FindYIndexNotLess(y, i_l, i_mid);
    }
    unsigned int FindYIndexNotLess(double y) { return FindYIndexNotLess(y, 0, NY - 1); }
    double Y0ForPositiveFAndNegativeY(double f, double alpha, double ymax, double fmax) {
        assert(f >= 0);
        assert(alpha < 7.97884560802865355E-1);
        // Recall that f(y) is downward sloping: 'left' and 'right' here refers to the ordering in f, i.e., the inverted direction of the ordering in y.
        double y, f_r = (Y[0] < ymax) ? fmax : F(0, alpha), f_l = 0;
        if (f >= f_r)
            y = (f - LN_TWO) / alpha;
        else {
            const unsigned int k = FindInterval(f, alpha, ymax, fmax, f_r, f_l);
            const double y_r = Y[k], y_l = Y[k + 1], dydf_r = 1 / Fp(k, alpha), dydf_l = 1 / Fp(k + 1u, alpha), y_i = G[k + 1], f_i = Fi(k + 1, alpha);
            const double r = InterpolationRationalCubic::ControlParameterToFitInteriorPoint(f_l, f_r, y_l, y_r, dydf_l, dydf_r, f_i, y_i);
            y = rational_cubic_interpolation(f, f_l, f_r, y_l, y_r, dydf_l, dydf_r, r);
        }
        return y;
    }
    double Y0(double f, double alpha) {
        double y;
        // Generic Maxima helper code used in parts of the below expansions:
        //   reset()$ kill(all)$ forget(facts())$ fpprec:75$ fpprintprec:18$ powerdisp:true$ ratprint:false$ load(revert)$ load(powers)$
        //   replace_derivatives(expr,f,x,d):=block([tmp,i],local(tmp,i),tmp:expr,for i:2 thru d do (tmp:subst(concat(f,i,p),'diff(f(x),x,i),tmp)),tmp:subst(fp,'diff(f(x),x),tmp),fullratsimp(subst(f,f(x),tmp)))$
        //   replace_derivatives_for_householder(expr,f,x,d):=block([tmp,i],local(tmp,i),tmp:expr,for i:2 thru d do (tmp:subst(concat(h,i)*fp,'diff(f(x),x,i),tmp)),tmp:subst(fp,'diff(f(x),x),tmp),fullratsimp(subst(-nu*fp,f(x),tmp)))$
        //   myhornerform(e,x):= block([f,s,k],local(f,s,k),f:expand(e),k:lmax(powers(f,x)),s:string(coeff(f,x,0)), for i:1 thru k do ( s:sconcat(s,"+",x,"*(",string(coeff(f,x,i))) ), for i:1 thru k do ( s:sconcat(s,")") ), eval_string(s) )$
        //   rathorner(expr,x):=block([tmp],local(tmp), tmp:coeff(taylor(denom(expr,x),x,0,1),x,0), (myhornerform(expand(num(expr)/tmp),x))/(myhornerform(expand(denom(expr)/tmp),x)))$
        //   filterbfloat(expr):=ssubst("(1-","(1.0E0-",ssubst("(1+","(1.0E0+",ssubst("(","(1.0E0*",ssubst("-","-1.0E0*",ssubst("+","+1.0E0*",ssubst("E","b",string(expr)))))))$
        //   /* nu = -f/f', h2=f''/f', h3=f'''/f', ... , hn=f^(n)/f' */
        //   HouseHolder(d):=block([f,x],local(f,d),rathorner(replace_derivatives_for_householder(d*(diff(1/f(x),x,d-1))/(diff(1/f(x),x,d)),f,x,d),nu))$
        if (f < 0 && alpha < 0.1388) { // alpha≈0.1388 <=> ymax≈-1.5
            // Here, y>0, and for increasing y, very rapidly, f(y) ≈ ln(2/π)/2-y²/2+α·y, i.e., essentially quadratic.
            // We approximate the inverse as a square root with offset, matching the behaviour near zero to fifth order.
            //   m:5$ r(nu):=nu*(1+nu*(a2+nu*a3))/(1+nu*(b1+nu*b2)); sqrt_:taylor(sqrt(1+2*r(nu))-1,nu,0,m); hh_:taylor(HouseHolder(5),nu,0,m);
            //   sols:solve([ coeff(sqrt_,nu,2)=coeff(hh_,nu,2), coeff(sqrt_,nu,3)=coeff(hh_,nu,3), coeff(sqrt_,nu,4)=coeff(hh_,nu,4), coeff(sqrt_,nu,5)=coeff(hh_,nu,5) ],[a2,a3,b1,b2])[1];
            // gives
            //   a2=-(-45*h2^2-45*h2^3+225*h2^4-135*h2^5+(60+90*h2-300*h2^2+150*h2^3)*h3+(40+20*h2)*h3^2+(-15+60*h2-45*h2^2-10*h3)*h4+(-6+6*h2)*h5)/(90*h2^2-180*h2^3+90*h2^4+(-120+180*h2-60*h2^2)*h3-40*h3^2+(-30+30*h2)*h4)
            //   a3=(270*h2^3-405*h2^4+135*h2^6+(-360*h2+630*h2^2-270*h2^4)*h3+(-180+180*h2^2)*h3^2-80*h3^3+(90-90*h2+60*h2*h3)*h4-15*h4^2+(18-18*h2^2+12*h3)*h5)/(540*h2^2-1080*h2^3+540*h2^4+(-720+1080*h2-360*h2^2)*h3-240*h3^2+(-180+180*h2)*h4)
            //   b1=-(-90*h2^3+180*h2^4-90*h2^5+(120*h2-210*h2^2+90*h2^3)*h3+(10+20*h2)*h3^2+(-15+45*h2-30*h2^2-5*h3)*h4+(-3+3*h2)*h5)/(45*h2^2-90*h2^3+45*h2^4+(-60+90*h2-30*h2^2)*h3-20*h3^2+(-15+15*h2)*h4)
            //   b2=(135*h2^4-270*h2^5+135*h2^6+(-120*h2^2+300*h2^3-180*h2^4)*h3+(-80+40*h2+40*h2^2)*h3^2-40*h3^3+(60*h2-90*h2^2+30*h2^3+(-20+40*h2)*h3)*h4-5*h4^2+(12*h2-12*h2^2+4*h3)*h5)/(180*h2^2-360*h2^3+180*h2^4+(-240+360*h2-120*h2^2)*h3-80*h3^2+(-60+60*h2)*h4)
            // which we then manually reorganize as below.
            const double f2p0 = -6.36619772367581343E-1, f3p0 = -2.1801361414499016E-1, f4p0 = 1.14770682054218857E-1, f5p0 = 4.43768846261782098E-3; // For reference: f6p0 = -1.5265853384165693E-1;
            const double fp0 = -7.97884560802865355E-1 + alpha, nu = f / fp0, h2 = f2p0 / fp0, h3 = f3p0 / fp0, h4 = f4p0 / fp0, h5 = f5p0 / fp0, h22 = h2 * h2, h23 = h22 * h2, h24 = h23 * h2, omh22 = (1 - h2) * (1 - h2), h32 = h3 * h3, h33 = h32 * h3, h42 = h4 * h4;
            const double a1 = 90 * h22 * omh22 - 60 * (2 - 3 * h2 + h22) * h3 - 40 * h32 + 30 * (h2 - 1) * h4;
            const double a2 = 45 * h22 * (1 + 3 * h2) * omh22 - 30 * (2 + 3 * h2 - 10 * h22 + 5 * h23) * h3 - 20 * (2 + h2) * h32 + (15 - 60 * h2 + 45 * h22 + 10 * h3) * h4 + 6 * (1 - h2) * h5;
            const double a3 = 135 * h23 * (2 + h2) * omh22 - 90 * (4 - 7 * h2 + 3 * h23) * h2 * h3 + 180 * (h22 - 1) * h32 - 80 * h33 + (90 - 90 * h2 + 60 * h2 * h3) * h4 - 15 * h42 + (18 - 18 * h22 + 12 * h3) * h5;
            const double b1 = 90 * h23 * omh22 - 30 * (4 - 7 * h2 + 3 * h22) * h2 * h3 - 10 * (1 + 2 * h2) * h32 + (15 - 45 * h2 + 30 * h22 + 5 * h3) * h4 + 3 * (1 - h2) * h5;
            const double b2 = 135 * h24 * omh22 - 60 * (2 - 5 * h2 + 3 * h22) * h22 * h3 - 40 * (2 - h2 - h22) * h32 - 40 * h33 + 10 * (6 * h2 - 9 * h22 + 3 * h23 - 2 * (1 - 2 * h2) * h3) * h4 - 5 * h42 + 4 * (3 * h2 - 3 * h22 + h3) * h5;
            y = SqrtOnePlusXMinusOne(2 * nu * (6 * a1 + nu * (6 * a2 + nu * a3)) / (6 * a1 + nu * (12 * b1 + 3 * nu * b2)));
        } else {
            if (alpha <= 0) { // Put option branch for f>=0 (the case f<0 is already handled above).
                y = Y0ForPositiveFAndNegativeY(f, alpha, -DBL_MAX, DBL_MAX);
            } else {
                // We know that at the maximum, fmax = ln(2·Φ(-ymax))+α·ymax = -ln(α/√(2/π))-ln(α)+α·ymax-ymax²/2 .
                // The latter expression becomes useful particularly when α ≈ √(2/π), since then we have y ≈ 0,
                // and the reformulation fmax = = -ln(α/√(2/π))-ln(α)+α·ymax-ymax²/2 then avoids subtractive
                // cancellation even without the need for any local expansions. In theory...
                // In practice, we found that fmax given by (-ln(α/√(2/π))-ln(α)+α·ymax-ymax²/2) can still be less than
                // fmax given by ln(2·Φ(-ymax))+α·ymax. For maximum consistency, we compute fmax from ymax exactly how
                // we compute f(y) in the subsequent iteration.
                double ymax = SQRT_TWO * InverseLnErfcx(log(SQRT_TWO_OVER_PI / alpha)), fmax = LnTwoCdf(-ymax) + alpha * ymax;
                if (fabs(f - fmax) <= DBL_EPSILON * fabs(fmax))
                    return ymax;
                if (f > fmax) {
                    fmax = LnTwoCdf(-ymax) + alpha * ymax;
                    if (fabs(f - fmax) <= DBL_EPSILON * fabs(fmax))
                        return ymax;
                    return NaN;
                }
                // The location of the maximum. fmax = f(ymax,alpha) = log(SqrtTwoOverPi/alpha)+alpha*ymax-0.5*ymax*ymax.
                // At the maximum, q = phi(-y)/Phi(-y) = α !
                if (alpha < 0.1388) { // 0 < α < 0.1388  ==>  ymax ≲ -1.5
                    // The case f<0 is already handled above
                    assert(f >= 0);
                    assert(ymax < -1);
                    double ymax_r, f_ymax_r, fp_ymax_r, f2p_ymax_r;
                    unsigned int i_ymax_r = FindYIndexNotLess(ymax);
                    assert(i_ymax_r + 1 < NY); // Because ymax ≲ -1.5.
                    const double y_shift = 0.5 * (Y[i_ymax_r + 1] - Y[i_ymax_r]);
                    i_ymax_r = FindYIndexNotLess(ymax + y_shift);
                    f_ymax_r = F(i_ymax_r, alpha);
                    if (f > f_ymax_r) {
                        // Transformed rational cubic interpolation of the inverse between the maximum (which is to the left of 0, i.e., ymax<0), and hl=zero (corresponding to ymax).
                        // Specifically, define:
                        //   w := y-ymax
                        //   h := sqrt(fmax-f)/r2
                        //   w = w(h)
                        //   yl = ymax
                        //   ->   fl = fmax,  wl = w(hl) = 0,       hl = 0,                 w'(hl)=1
                        //   yr = - ⌊ |ymax| ⌋
                        //   ->   fr = f(yr), wr = w(hr) = yr-ymax, hr=sqrt(fmax-f(yr))/r2, w'(hr)=-2·r2·sqrt(fmax-f(yr))/f'(yr),   w''(hr)=-2·r2²·[1+2·(fmax-f(yr))·f''(yr)/f'(yr)²]/f'(yr)
                        // Rational cubic interpolation of w(h) over h \in [hl,hr]
                        ymax_r = Y[i_ymax_r];
                        fp_ymax_r = Fp(i_ymax_r, alpha);
                        f2p_ymax_r = E[i_ymax_r];
                        const double q = alpha, c2 = q * (ymax - q) / 2, r2 = sqrt(fabs(c2)), h = sqrt(fmax - f) / r2;
                        const double sqrt_df = sqrt(fmax - f_ymax_r), hr = sqrt_df / r2, wr = ymax_r - ymax, sr = -2 * r2 * sqrt_df / fp_ymax_r;
                        const double r = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(0, hr, 0, wr, 1, sr, -2 * r2 * r2 / fp_ymax_r * (1 + 2 * (fmax - f_ymax_r) * (f2p_ymax_r / fp_ymax_r) / fp_ymax_r), true);
                        y = ymax + rational_cubic_interpolation(h, 0, hr, 0, wr, 1, sr, r);
                    } else
                        y = Y0ForPositiveFAndNegativeY(f, alpha, ymax, fmax);
                } else {
                    // Define:
                    //   w := y-ymax
                    //   h := sqrt(fmax-f)/r2
                    // Expand h(w) as a Taylor series around y=ymax, i.e., w=0. Since f'(ymax)=0, this expansion is linear in w in leading order.
                    // Using the Lagrange inversion theorem, invert this expansion as w=w(h), now linear in h in leading order. Tidy up. In short:
                    //   tmp:1+w*(r3+w*(r4+w*r5)); tmp:taylor(w*r2*sqrt(1+w*(r3+w*(r4+w*r5))),w,0,3); tmp:revert(tmp,w); tmp:ev(tmp,w=h*r2); horner(tmp,h); tmp:pade(taylor(tmp,h,0,3),2,1)[1]; rathorner(tmp,h);
                    // First, however, we need the coefficients of the chosen form h(w) ≈ w*r2*sqrt(1+w*(r3+w*(r4+w*r5))).
                    // These are derived from the the Taylor series of f(y) = fmax+y*y*(c2+y*(c3+y*(c4+y*(c5...)))) with ci = f^(i)(ymax)/i! .
                    double q = alpha, c2 = q * (ymax - q) / 2, r2 = sqrt(fabs(c2)), h = sqrt(fmax - f) / r2;
                    double ymax2 = ymax * ymax, r3 = (q * (1 + q * (-2 * q + 3 * ymax) - ymax2)) / 6 / c2, r4 = (q * (q * (4 + q * (-6 * q + 12 * ymax) - 7 * ymax2) + ymax * (-3 + ymax2))) / 24 / c2; // For reference: r5 = (q * (-3 + ymax2 * (6 - ymax2) + q * (q * (20 + q * (-24 * q + 60 * ymax) - 50 * ymax2) + ymax * (-25 + 15 * ymax2)))) / 120 / c2;
                    double R = (h * (1 + h * ((3 * r3) / 4 - r4 / r3))) / (1 + h * ((5 * r3) / 4 - r4 / r3));
                    if (ymax < 0 && f > -10 * fmax) { // Adjust the expansion out of the maximum such that it goes exactly through the origin at (y=0,f=0), but only if the maximum is to the left of.
                        assert(fmax >= 0);
                        const double h0 = sqrt(fmax - 0) / r2, R0 = (h0 * (1 + h0 * ((3 * r3) / 4 - r4 / r3))) / (1 + h0 * ((5 * r3) / 4 - r4 / r3));
                        y = (f * (fmax * (ymax + R0) + f * (ymax + R)) + fmax * fmax * (R - R0)) / (fmax * fmax + f * f);
                    } else
                        y = ymax + R;
                }
            }
        }
        // Improvement of initial guess based on asymptotic approximation for y -> +infty where f(y) ~ ln(2/pi)/2+alpha·y-y²/2-ln(y) + O(1/y²)
        // This is applied only when f<0 and y>y_threshold for some moderately large y_threshold. We use y_threshold=5.
        // Then, we solve the quadratic expansion around the current initial guess y. Note that ln(2/pi) = -0.45158270528945486
        // Maxima code:   tmp:ev(log(2/%pi)/2-f+alpha*y-y^2/2-log(y),y=yp*(1+t));taylor(tmp,t,0,2);
        // gives:         -(2*log(yp)-log(2/%pi)+2*f-2*alpha*yp+yp^2)/2+(-1+alpha*yp-yp^2)*t-((-1+yp^2)*t^2)/2+...
        // This "pre-polishing" enables the method to converge for very large values of y. Without this, when, say, alpha=1 and y=300, the above initial guess does not enter the domain of attraction of the subsequent Householder iteration.
        if (f < 0 && y > 5) {
            const double a = (1 - y * y) / 2, b = y * (alpha - y) - 1, c = y * (alpha - y / 2) - log(y) - 0.45158270528945486 / 2 - f, p = b / a;
            y *= 1 - 0.5 * p + sqrt(0.25 * p * p - c / a);
        }
        return y;
    }
    double HouseHolder6Improvement(double y, double f, double alpha) {
        const double ln2cdf = LnTwoCdf(-y), q = exp(-0.5 * y * y - ln2cdf) * SQRT_TWO_OVER_PI, fy = ln2cdf + alpha * y, f1 = alpha - q, f2 = q * (y - q), y2 = y * y, f3 = (q * (1 + q * (-2 * q + 3 * y) - y2)), f4 = (q * (q * (4 + q * (-6 * q + 12 * y) - 7 * y2) + y * (-3 + y2))), f5 = (q * (-3 + y2 * (6 - y2) + q * (q * (20 + q * (-24 * q + 60 * y) - 50 * y2) + y * (-25 + 15 * y2)))), f6 = q * (y * (15 + y2 * (-10 + y2)) + q * (-28 + y2 * (101 - 31 * y2) + q * (q * (120 + q * (-120 * q + 360 * y) - 390 * y2) + y * (-210 + 180 * y2))));
        return HouseHolder6Step(-(fy - f) / f1, f2 / f1, f3 / f1, f4 / f1, f5 / f1, f6 / f1);
    }
}

// f(y) = ln(2·Φ(-y)) + α·y  as defined in equation (1.14)
EXPORT_EXTERN_C double ModifiedDeltaWithPremiumFunction(double y, double alpha) { return LnTwoCdf(-y) + alpha * y; }

DECLARE_XL_FUNCTION( ModifiedDeltaWithPremiumFunction, 0 , "BBB" , "y,alpha", "returns the modified delta-with-premium function f(y;alpha) as defined in formula (1.14).", "the transformed moneyness.", "the delta-signed standard deviation")

EXPORT_EXTERN_C double InverseModifiedDeltaWithPremiumFunction(double f, double alpha) {
    const double y = Y0(f, alpha);
    // ONE HouseHolder(6) improvement suffices for all inputs.
    return y + HouseHolder6Improvement(y, f, alpha);
}

DECLARE_XL_FUNCTION( InverseModifiedDeltaWithPremiumFunction, 0 , "BBB" , "f,alpha", "returns the inverse modified delta-with-premium function y(f;alpha), the inverse of f(y;alpha) as defined in formula (1.14).", "the modfied delta-with-premium function value.", "the delta-signed standard deviation.")

double AttainableAccuracyFactor(double y, double alpha){ // max( |f(y)/(y·f'(y))| , 1 )
    return std::max( fabs(ModifiedDeltaWithPremiumFunction(y,alpha)/(y*(alpha-norm_pdf(y)/norm_cdf(-y)))) , 1.0 );
}
