//+------------------------------------------------------------------+
//|                                                         ndtr.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include "error.mqh"
#include "const.mqh"
#include "polevl.mqh"
namespace special {
namespace cephes {

    namespace detail {

         double ndtr_P[] = {2.46196981473530512524E-10, 5.64189564831068821977E-1, 7.46321056442269912687E0,
                                     4.86371970985681366614E1,   1.96520832956077098242E2,  5.26445194995477358631E2,
                                     9.34528527171957607540E2,   1.02755188689515710272E3,  5.57535335369399327526E2};

         double ndtr_Q[] = {
            /* 1.00000000000000000000E0, */
            1.32281951154744992508E1, 8.67072140885989742329E1, 3.54937778887819891062E2, 9.75708501743205489753E2,
            1.82390916687909736289E3, 2.24633760818710981792E3, 1.65666309194161350182E3, 5.57535340817727675546E2};

         double ndtr_R[] = {5.64189583547755073984E-1, 1.27536670759978104416E0, 5.01905042251180477414E0,
                                     6.16021097993053585195E0,  7.40974269950448939160E0, 2.97886665372100240670E0};

         double ndtr_S[] = {
            /* 1.00000000000000000000E0, */
            2.26052863220117276590E0, 9.39603524938001434673E0, 1.20489539808096656605E1,
            1.70814450747565897222E1, 9.60896809063285878198E0, 3.36907645100081516050E0};

         double ndtr_T[] = {9.60497373987051638749E0, 9.00260197203842689217E1, 2.23200534594684319226E3,
                                     7.00332514112805075473E3, 5.55923013010394962768E4};

         double ndtr_U[] = {
            /* 1.00000000000000000000E0, */
            3.35617141647503099647E1, 5.21357949780152679795E2, 4.59432382970980127987E3, 2.26290000613890934246E4,
            4.92673942608635921086E4};

         double ndtri_UTHRESH = 37.519379347;

    } // namespace detail


     inline double erfc(double a) {
        double p, q, x, y, z;

        if (isnan(a)) {
            set_error("erfc", SF_ERROR_DOMAIN, NULL);
            return quiet_NaN();
        }

        if (a < 0.0) {
            x = -a;
        } else {
            x = a;
        }

        if (x < 1.0) {
            return 1.0 - erf(a);
        }

        z = -a * a;

        if (z < -detail::MAXLOG) {
            set_error("erfc", SF_ERROR_UNDERFLOW, NULL);
        if (a < 0) {
            return 2.0;
        } else {
            return 0.0;
        }
        }

        z = exp(z);

        if (x < 8.0) {
            p = polevl(x, detail::ndtr_P, 8);
            q = p1evl(x, detail::ndtr_Q, 8);
        } else {
            p = polevl(x, detail::ndtr_R, 5);
            q = p1evl(x, detail::ndtr_S, 6);
        }
        y = (z * p) / q;

        if (a < 0) {
            y = 2.0 - y;
        }

        if (y != 0.0) {
            return y;
        }

        set_error("erfc", SF_ERROR_UNDERFLOW, NULL);
        if (a < 0) {
            return 2.0;
        } else {
            return 0.0;
        }
    }

     inline double erf(double x) {
        double y, z;

        if (isnan(x)) {
            set_error("erf", SF_ERROR_DOMAIN, NULL);
            return quiet_NaN();
        }

        if (x < 0.0) {
            return -erf(-x);
        }

        if (abs(x) > 1.0) {
            return (1.0 - erfc(x));
        }
        z = x * x;

        y = x * polevl(z, detail::ndtr_T, 4) / p1evl(z, detail::ndtr_U, 5);
        return y;
    }

     inline double ndtr(double a) {
        double x, y, z;

        if (isnan(a)) {
            set_error("ndtr", SF_ERROR_DOMAIN, NULL);
            return quiet_NaN();
        }

        x = a * M_SQRT1_2;
        z = abs(x);

        if (z < M_SQRT1_2) {
            y = 0.5 + 0.5 * erf(x);
        } else {
            y = 0.5 * erfc(z);
            if (x > 0) {
                y = 1.0 - y;
            }
        }

        return y;
    }

} // namespace cephes
} // namespace special