/**
 * Error Functions and Normal Distribution.
 *
 * License: $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0).
 * Copyright: Based on the CEPHES math library, which is
 *            Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
 * Authors:   Stephen L. Moshier, ported to D by Don Clugston and David Nadlinger
 */
/**
 * Macros:
 *  NAN = $(RED NAN)
 *  SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
 *  GAMMA =  &#915;
 *  INTEGRAL = &#8747;
 *  INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
 *  POWER = $1<sup>$2</sup>
 *  BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
 *  CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
 *  TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
 *      <caption>Special Values</caption>
 *      $0</table>
 *  SVH = $(TR $(TH $1) $(TH $2))
 *  SV  = $(TR $(TD $1) $(TD $2))
 */
module std.internal.math.errorfunction;
import std.math;

pure:
nothrow:
@safe:
@nogc:

private {
immutable real EXP_2  = 0.135335283236612691893999494972484403L; /* exp(-2) */
enum real SQRT2PI = 2.50662827463100050241576528481104525L; // sqrt(2pi)


enum real MAXLOG = 0x1.62e42fefa39ef358p+13L;  // log(real.max)
enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min*real.epsilon) = log(smallest denormal)
}

T rationalPoly(T)(T x, const(T) [] numerator, const(T) [] denominator) pure nothrow
{
    return poly(x, numerator)/poly(x, denominator);
}


private {

/* erfc(x) = exp(-x^2) P(1/x)/Q(1/x)
   1/8 <= 1/x <= 1
   Peak relative error 5.8e-21  */
immutable real[10] P = [ -0x1.30dfa809b3cc6676p-17, 0x1.38637cd0913c0288p+18,
   0x1.2f015e047b4476bp+22, 0x1.24726f46aa9ab08p+25, 0x1.64b13c6395dc9c26p+27,
   0x1.294c93046ad55b5p+29, 0x1.5962a82f92576dap+30, 0x1.11a709299faba04ap+31,
   0x1.11028065b087be46p+31, 0x1.0d8ef40735b097ep+30
];

immutable real[11] Q = [ 0x1.14d8e2a72dec49f4p+19, 0x1.0c880ff467626e1p+23,
   0x1.04417ef060b58996p+26, 0x1.404e61ba86df4ebap+28, 0x1.0f81887bc82b873ap+30,
   0x1.4552a5e39fb49322p+31, 0x1.11779a0ceb2a01cep+32, 0x1.3544dd691b5b1d5cp+32,
   0x1.a91781f12251f02ep+31, 0x1.0d8ef3da605a1c86p+30, 1.0
];

// For 128 bit quadruple-precision floats, we use a higher-precision implementation
// with more polynomial segments.
enum isIEEEQuadruple = floatTraits!real.realFormat == RealFormat.ieeeQuadruple;
static if (isIEEEQuadruple)
{
    // erfc(x + 0.25) = erfc(0.25) + x R(x)
    // 0 <= x < 0.125
    // Peak relative error 1.4e-35
    immutable real[9] RNr13 = [
        -2.353707097641280550282633036456457014829E3L,
        3.871159656228743599994116143079870279866E2L,
        -3.888105134258266192210485617504098426679E2L,
        -2.129998539120061668038806696199343094971E1L,
        -8.125462263594034672468446317145384108734E1L,
        8.151549093983505810118308635926270319660E0L,
        -5.033362032729207310462422357772568553670E0L,
        -4.253956621135136090295893547735851168471E-2L,
        -8.098602878463854789780108161581050357814E-2L
    ];
    immutable real[9] RDr13 = [
        2.220448796306693503549505450626652881752E3L,
        1.899133258779578688791041599040951431383E2L,
        1.061906712284961110196427571557149268454E3L,
        7.497086072306967965180978101974566760042E1L,
        2.146796115662672795876463568170441327274E2L,
        1.120156008362573736664338015952284925592E1L,
        2.211014952075052616409845051695042741074E1L,
        6.469655675326150785692908453094054988938E-1L,
        1.0
    ];

    // erfc(0.25) = C13a + C13b to extra precision.
    immutable real C13a = 0.723663330078125L;
    immutable real C13b = 1.0279753638067014931732235184287934646022E-5L;

    // erfc(x + 0.375) = erfc(0.375) + x R(x)
    // 0 <= x < 0.125
    // Peak relative error 1.2e-35
    immutable real[9] RNr14 = [
        -2.446164016404426277577283038988918202456E3L,
        6.718753324496563913392217011618096698140E2L,
        -4.581631138049836157425391886957389240794E2L,
        -2.382844088987092233033215402335026078208E1L,
        -7.119237852400600507927038680970936336458E1L,
        1.313609646108420136332418282286454287146E1L,
        -6.188608702082264389155862490056401365834E0L,
        -2.787116601106678287277373011101132659279E-2L,
        -2.230395570574153963203348263549700967918E-2L
    ];
    immutable real[9] RDr14 = [
        2.495187439241869732696223349840963702875E3L,
        2.503549449872925580011284635695738412162E2L,
        1.159033560988895481698051531263861842461E3L,
        9.493751466542304491261487998684383688622E1L,
        2.276214929562354328261422263078480321204E2L,
        1.367697521219069280358984081407807931847E1L,
        2.276988395995528495055594829206582732682E1L,
        7.647745753648996559837591812375456641163E-1L,
        1.0
    ];

    // erfc(0.375) = C14a + C14b to extra precision.
    immutable real C14a = 0.5958709716796875L;
    immutable real C14b = 1.2118885490201676174914080878232469565953E-5L;

    // erfc(x + 0.5) = erfc(0.5) + x R(x)
    // 0 <= x < 0.125
    // Peak relative error 4.7e-36
    immutable real[9] RNr15 = [
        -2.624212418011181487924855581955853461925E3L,
        8.473828904647825181073831556439301342756E2L,
        -5.286207458628380765099405359607331669027E2L,
        -3.895781234155315729088407259045269652318E1L,
        -6.200857908065163618041240848728398496256E1L,
        1.469324610346924001393137895116129204737E1L,
        -6.961356525370658572800674953305625578903E0L,
        5.145724386641163809595512876629030548495E-3L,
        1.990253655948179713415957791776180406812E-2L
    ];
    immutable real[9] RDr15 = [
        2.986190760847974943034021764693341524962E3L,
        5.288262758961073066335410218650047725985E2L,
        1.363649178071006978355113026427856008978E3L,
        1.921707975649915894241864988942255320833E2L,
        2.588651100651029023069013885900085533226E2L,
        2.628752920321455606558942309396855629459E1L,
        2.455649035885114308978333741080991380610E1L,
        1.378826653595128464383127836412100939126E0L,
        1.0
    ];
    // erfc(0.5) = C15a + C15b to extra precision.
    immutable real C15a = 0.4794921875L;
    immutable real C15b = 7.9346869534623172533461080354712635484242E-6L;

    // erfc(x + 0.625) = erfc(0.625) + x R(x)
    // 0 <= x < 0.125
    // Peak relative error 5.1e-36
    immutable real[9] RNr16 = [
        -2.347887943200680563784690094002722906820E3L,
        8.008590660692105004780722726421020136482E2L,
        -5.257363310384119728760181252132311447963E2L,
        -4.471737717857801230450290232600243795637E1L,
        -4.849540386452573306708795324759300320304E1L,
        1.140885264677134679275986782978655952843E1L,
        -6.731591085460269447926746876983786152300E0L,
        1.370831653033047440345050025876085121231E-1L,
        2.022958279982138755020825717073966576670E-2L,
    ];
    immutable real[9] RDr16 = [
        3.075166170024837215399323264868308087281E3L,
        8.730468942160798031608053127270430036627E2L,
        1.458472799166340479742581949088453244767E3L,
        3.230423687568019709453130785873540386217E2L,
        2.804009872719893612081109617983169474655E2L,
        4.465334221323222943418085830026979293091E1L,
        2.612723259683205928103787842214809134746E1L,
        2.341526751185244109722204018543276124997E0L,
        1.0
    ];
    // erfc(0.625) = C16a + C16b to extra precision.
    immutable real C16a = 0.3767547607421875L;
    immutable real C16b = 4.3570693945275513594941232097252997287766E-6L;

    // erfc(x + 0.75) = erfc(0.75) + x R(x)
    // 0 <= x < 0.125
    // Peak relative error 1.7e-35
    immutable real[9] RNr17 = [
        -1.767068734220277728233364375724380366826E3L,
        6.693746645665242832426891888805363898707E2L,
        -4.746224241837275958126060307406616817753E2L,
        -2.274160637728782675145666064841883803196E1L,
        -3.541232266140939050094370552538987982637E1L,
        6.988950514747052676394491563585179503865E0L,
        -5.807687216836540830881352383529281215100E0L,
        3.631915988567346438830283503729569443642E-1L,
        -1.488945487149634820537348176770282391202E-2L
    ];
    immutable real[9] RDr17 = [
        2.748457523498150741964464942246913394647E3L,
        1.020213390713477686776037331757871252652E3L,
        1.388857635935432621972601695296561952738E3L,
        3.903363681143817750895999579637315491087E2L,
        2.784568344378139499217928969529219886578E2L,
        5.555800830216764702779238020065345401144E1L,
        2.646215470959050279430447295801291168941E1L,
        2.984905282103517497081766758550112011265E0L,
        1.0
    ];
    // erfc(0.75) = C17a + C17b to extra precision.
    immutable real C17a = 0.2888336181640625L;
    immutable real C17b = 1.0748182422368401062165408589222625794046E-5L;


    // erfc(x + 0.875) = erfc(0.875) + x R(x)
    // 0 <= x < 0.125
    // Peak relative error 2.2e-35
    immutable real[9] RNr18 = [
        -1.342044899087593397419622771847219619588E3L,
        6.127221294229172997509252330961641850598E2L,
        -4.519821356522291185621206350470820610727E2L,
        1.223275177825128732497510264197915160235E1L,
        -2.730789571382971355625020710543532867692E1L,
        4.045181204921538886880171727755445395862E0L,
        -4.925146477876592723401384464691452700539E0L,
        5.933878036611279244654299924101068088582E-1L,
        -5.557645435858916025452563379795159124753E-2L
    ];
    immutable real[9] RDr18 = [
        2.557518000661700588758505116291983092951E3L,
        1.070171433382888994954602511991940418588E3L,
        1.344842834423493081054489613250688918709E3L,
        4.161144478449381901208660598266288188426E2L,
        2.763670252219855198052378138756906980422E2L,
        5.998153487868943708236273854747564557632E1L,
        2.657695108438628847733050476209037025318E1L,
        3.252140524394421868923289114410336976512E0L,
        1.0
    ];

    // erfc(0.875) = C18a + C18b to extra precision.
    immutable real C18a = 0.215911865234375L;
    immutable real C18b = 1.3073705765341685464282101150637224028267E-5L;

    // erfc(x + 1.0) = erfc(1.0) + x R(x)
    // 0 <= x < 0.125
    // Peak relative error 1.6e-35
    immutable real[9] RNr19 = [
        -1.139180936454157193495882956565663294826E3L,
        6.134903129086899737514712477207945973616E2L,
        -4.628909024715329562325555164720732868263E2L,
        4.165702387210732352564932347500364010833E1L,
        -2.286979913515229747204101330405771801610E1L,
        1.870695256449872743066783202326943667722E0L,
        -4.177486601273105752879868187237000032364E0L,
        7.533980372789646140112424811291782526263E-1L,
        -8.629945436917752003058064731308767664446E-2L
    ];
    immutable real[9] RDr19 = [
        2.744303447981132701432716278363418643778E3L,
        1.266396359526187065222528050591302171471E3L,
        1.466739461422073351497972255511919814273E3L,
        4.868710570759693955597496520298058147162E2L,
        2.993694301559756046478189634131722579643E2L,
        6.868976819510254139741559102693828237440E1L,
        2.801505816247677193480190483913753613630E1L,
        3.604439909194350263552750347742663954481E0L,
        1.0
    ];

    // erfc(1.0) = C19a + C19b to extra precision.
    immutable real C19a = 0.15728759765625L;
    immutable real C19b = 1.1609394035130658779364917390740703933002E-5L;

    // erfc(x + 1.125) = erfc(1.125) + x R(x)
    // 0 <= x < 0.125
    // Peak relative error 3.6e-36
    immutable real[9] RNr20 = [
        -9.652706916457973956366721379612508047640E2L,
        5.577066396050932776683469951773643880634E2L,
        -4.406335508848496713572223098693575485978E2L,
        5.202893466490242733570232680736966655434E1L,
        -1.931311847665757913322495948705563937159E1L,
        -9.364318268748287664267341457164918090611E-2L,
        -3.306390351286352764891355375882586201069E0L,
        7.573806045289044647727613003096916516475E-1L,
        -9.611744011489092894027478899545635991213E-2L
    ];
    immutable real[9] RDr20 = [
        3.032829629520142564106649167182428189014E3L,
        1.659648470721967719961167083684972196891E3L,
        1.703545128657284619402511356932569292535E3L,
        6.393465677731598872500200253155257708763E2L,
        3.489131397281030947405287112726059221934E2L,
        8.848641738570783406484348434387611713070E1L,
        3.132269062552392974833215844236160958502E1L,
        4.430131663290563523933419966185230513168E0L,
        1.0
    ];

    // erfc(1.125) = C20a + C20b to extra precision.
    immutable real C20a = 0.111602783203125L;
    immutable real C20b = 8.9850951672359304215530728365232161564636E-6L;

    // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
    // 7/8 <= 1/x < 1
    // Peak relative error 1.4e-35
    immutable real[10] RNr8 = [
        3.587451489255356250759834295199296936784E1L,
        5.406249749087340431871378009874875889602E2L,
        2.931301290625250886238822286506381194157E3L,
        7.359254185241795584113047248898753470923E3L,
        9.201031849810636104112101947312492532314E3L,
        5.749697096193191467751650366613289284777E3L,
        1.710415234419860825710780802678697889231E3L,
        2.150753982543378580859546706243022719599E2L,
        8.740953582272147335100537849981160931197E0L,
        4.876422978828717219629814794707963640913E-2L
    ];
    immutable real[10] RDr8 = [
        6.358593134096908350929496535931630140282E1L,
        9.900253816552450073757174323424051765523E2L,
        5.642928777856801020545245437089490805186E3L,
        1.524195375199570868195152698617273739609E4L,
        2.113829644500006749947332935305800887345E4L,
        1.526438562626465706267943737310282977138E4L,
        5.561370922149241457131421914140039411782E3L,
        9.394035530179705051609070428036834496942E2L,
        6.147019596150394577984175188032707343615E1L,
        1.0L
    ];

    // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
    // 3/4 <= 1/x < 7/8
    // Peak relative error 1.7e-36
    immutable real[10] RNr7 = [
        1.293515364043117601705535485785956592493E2L,
        2.474534371269189867053251150063459055230E3L,
        1.756537563959875738809491329250457486510E4L,
        5.977479535376639763773153344676726091607E4L,
        1.054313816139671870123172936972055385049E5L,
        9.754699773487726957401038094714603033904E4L,
        4.579131242577671038339922925213209214880E4L,
        1.000710322164510887997115157797717324370E4L,
        8.496863250712471449526805271633794700452E2L,
        1.797349831386892396933210199236530557333E1L
    ];
    immutable real[11] RDr7 = [
        2.292696320307033494820399866075534515002E2L,
        4.500632608295626968062258401895610053116E3L,
        3.321218723485498111535866988511716659339E4L,
        1.196084512221845156596781258490840961462E5L,
        2.287033883912529843927983406878910939930E5L,
        2.370223495794642027268482075021298394425E5L,
        1.305173734022437154610938308944995159199E5L,
        3.589386258485887630236490009835928559621E4L,
        4.339996864041074149726360516336773136101E3L,
        1.753135522665469574605384979152863899099E2L,
        1.0L
    ];

    // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
    // 5/8 <= 1/x < 3/4
    // Peak relative error 1.6e-35
    immutable real[10] RNr6 = [
        1.423313561367394923305025174137639124533E1L,
        3.244462503273003837514629113846075327206E2L,
        2.784937282403293364911673341412846781934E3L,
        1.163060685810874867196849890286455473382E4L,
        2.554141394931962276102205517358731053756E4L,
        2.982733782500729530503336931258698708782E4L,
        1.789683564523810605328169719436374742840E4L,
        5.056032142227470121262177112822018882754E3L,
        5.605349942234782054561269306895707034586E2L,
        1.561652599080729507274832243665726064881E1L
    ];
    immutable real[11] RDr6 = [
        2.522757606611479946069351519410222913326E1L,
        5.876797910931896554014229647006604017806E2L,
        5.211092128250480712011248211246144751074E3L,
        2.282679910855404599271496827409168580797E4L,
        5.371245819205596609986320599133109262447E4L,
        6.926186102106400355114925675028888924445E4L,
        4.794366033363621432575096487724913414473E4L,
        1.673190682734065914573814938835674963896E4L,
        2.589544846151313120096957014256536236242E3L,
        1.349438432583208276883323156200117027433E2L,
        1.0L
    ];

    // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
    // 1/2 <= 1/x < 5/8
    // Peak relative error 4.3e-36
    immutable real[11] RNr5 = [
        6.743447478279267339131211137241149796763E-2L,
        2.031339015932962998168472743355874796350E0L,
        2.369234815713196480221800940201367484379E1L,
        1.387819614016107433603101545594790875922E2L,
        4.435600256936515839937720907171966121786E2L,
        7.881577949936817507981170892417739733047E2L,
        7.615749099291545976179905281851765734680E2L,
        3.752484528663442467089606663006771157777E2L,
        8.279644286027145214308303292537009564726E1L,
        6.201462983413738162709722770960040042647E0L,
        6.649631608720062333043506249503378282697E-2L
    ];
    immutable real[11] RDr5 = [
        1.195244945161889822018178270706903972343E-1L,
        3.660216908153253021384862427197665991311E0L,
        4.373405883243078019655721779021995159854E1L,
        2.653305963056235008916733402786877121865E2L,
        8.921329790491152046318422124415895506335E2L,
        1.705552231555600759729260640562363304312E3L,
        1.832711109606893446763174603477244625325E3L,
        1.056823953275835649973998168744261083316E3L,
        2.975561981792909722126456997074344895584E2L,
        3.393149095158232521894537008472203487436E1L,
        1.0L
    ];

    // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
    // 3/8 <= 1/x < 1/2
    // Peak relative error 1.8e-36
    immutable real[11] RNr4 = [
        3.558685919236420073872459554885612994007E-2L,
        1.460223642496950651561817195253277924528E0L,
        2.379856746555189546876720308841066577268E1L,
        2.005205521246554860334064698817220160117E2L,
        9.533374928664989955629120027419490517596E2L,
        2.623024576994438336130421711314560425373E3L,
        4.126446434603735586340585027628851620886E3L,
        3.540675861596687801829655387867654030013E3L,
        1.506037084891064572653273761987617394259E3L,
        2.630715699182706745867272452228891752353E2L,
        1.202476629652900619635409242749750364878E1L
    ];
    immutable real[12] RDr4 = [
        6.307606561714590590399683184410336583739E-2L,
        2.619717051134271249293056836082721776665E0L,
        4.344441402681725017630451522968410844608E1L,
        3.752891116408399440953195184301023399176E2L,
        1.849305988804654653921972804388006355502E3L,
        5.358505261991675891835885654499883449403E3L,
        9.091890995405251314631428721090705475825E3L,
        8.731418313949291797856351745278287516416E3L,
        4.420211285043270337492325400764271868740E3L,
        1.031487363021856106882306509107923200832E3L,
        8.387036084846046121805145056040429461783E1L,
        1.0L
    ];

    // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
    // 1/4 <= 1/x < 3/8
    // Peak relative error 8.1e-37
    immutable real[12] RNr3 = [
        4.584481542956275354582319313040418316755E-5L,
        2.674923158288848442110883948437930349128E-3L,
        6.344232532055212248017211243740466847311E-2L,
        7.985145965992002744933550450451513513963E-1L,
        5.845078061888281665064746347663123946270E0L,
        2.566625318816866587482759497608029522596E1L,
        6.736225182343446605268837827950856640948E1L,
        1.021796460139598089409347761712730512053E2L,
        8.344336615515430530929955615400706619764E1L,
        3.207749011528249356283897356277376306967E1L,
        4.386185123863412086856423971695142026036E0L,
        8.971576448581208351826868348023528863856E-2L
    ];
    immutable real[12] RDr3 = [
        8.125781965218112303281657065320409661370E-5L,
        4.781806762611504685247817818428945295520E-3L,
        1.147785538413798317790357996845767614561E-1L,
        1.469285552007088106614218996464752307606E0L,
        1.101712261349880339221039938999124077650E1L,
        5.008507527095093413713171655268276861071E1L,
        1.383058691613468970486425146336829447433E2L,
        2.264114250278912520501010108736339599752E2L,
        2.081377197698598680576330179979996940039E2L,
        9.724438129690013609440151781601781137944E1L,
        1.907905050771832372089975877589291760121E1L,
        1.0L
    ];

    // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
    // 1/8 <= 1/x < 1/4
    // Peak relative error 1.5e-36
    immutable real[12] RNr2 = [
        6.928615158005256885698045840589513728399E-7L,
        5.616245938942075826026382337922413007879E-5L,
        1.871624980715261794832358438894219696113E-3L,
        3.349922063795792371642023765253747563009E-2L,
        3.531865233974349943956345502463135695834E-1L,
        2.264714157625072773976468825160906342360E0L,
        8.810720294489253776747319730638214883026E0L,
        2.014056685571655833019183248931442888437E1L,
        2.524586947657190747039554310814128743320E1L,
        1.520656940937208886246188940244581671609E1L,
        3.334145500790963675035841482334493680498E0L,
        1.122108380007109245896534245151140632457E-1L
    ];
    immutable real[12] RDr2 = [
        1.228065061824874795984937092427781089256E-6L,
        1.001593999520159167559129042893802235969E-4L,
        3.366527555699367241421450749821030974446E-3L,
        6.098626947195865254152265585991861150369E-2L,
        6.541547922508613985813189387198804660235E-1L,
        4.301130233305371976727117480925676583204E0L,
        1.737155892350891711527711121692994762909E1L,
        4.206892112110558214680649401236873828801E1L,
        5.787487996025016843403524261574779631219E1L,
        4.094047148590822715163181507813774861621E1L,
        1.230603930056944875836549716515643997094E1L,
        1.0L
    ];

    // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
    // 1/128 <= 1/x < 1/8
    // Peak relative error 2.2e-36
    immutable real[10] RNr1 = [
        1.293111801138199795250035229383033322539E-6L,
        9.785224720880980456171438759161161816706E-5L,
        2.932474396233212166056331430621176065943E-3L,
        4.496429595719847083917435337780697436921E-2L,
        3.805989855927720844877478869846718877846E-1L,
        1.789745532222426292126781724570152590071E0L,
        4.465737379634389318903237306594171764628E0L,
        5.268535822258082278401240171488850433767E0L,
        2.258357766807433839494276681092713991651E0L,
        1.504459334078750002966538036652860809497E-1L
    ];
    immutable real[10] RDr1 = [
        2.291980991578770070179177302906728371406E-6L,
        1.745845828808028552029674694534934620384E-4L,
        5.283248841982102317072923869576785278019E-3L,
        8.221212297078141470944454807434634848018E-2L,
        7.120500674861902950423510939060230945621E-1L,
        3.475435367560809622183983439133664598155E0L,
        9.243253391989233533874386043611304387113E0L,
        1.227894792475280941511758877318903197188E1L,
        6.789361410398841316638617624392719077724E0L,
        1.0L
    ];

    // erf(z+1)  = erfConst + P(z)/Q(z)
    // -.125 <= z <= 0
    // Peak relative error 7.3e-36
    immutable real erfConst = 0.845062911510467529296875L;
    immutable real[9] TN2 = [
        -4.088889697077485301010486931817357000235E1L,
        7.157046430681808553842307502826960051036E3L,
        -2.191561912574409865550015485451373731780E3L,
        2.180174916555316874988981177654057337219E3L,
        2.848578658049670668231333682379720943455E2L,
        1.630362490952512836762810462174798925274E2L,
        6.317712353961866974143739396865293596895E0L,
        2.450441034183492434655586496522857578066E1L,
        5.127662277706787664956025545897050896203E-1L
    ];
    immutable real[10] TD2 = [
        1.731026445926834008273768924015161048885E4L,
        1.209682239007990370796112604286048173750E4L,
        1.160950290217993641320602282462976163857E4L,
        5.394294645127126577825507169061355698157E3L,
        2.791239340533632669442158497532521776093E3L,
        8.989365571337319032943005387378993827684E2L,
        2.974016493766349409725385710897298069677E2L,
        6.148192754590376378740261072533527271947E1L,
        1.178502892490738445655468927408440847480E1L,
        1.0L
    ];

    // erf(x)  = x + x P(x^2)/Q(x^2)
    // 0 <= x <= 7/8
    // Peak relative error 1.8e-35
    immutable real[9] TN1 = [
        -3.858252324254637124543172907442106422373E10L,
        9.580319248590464682316366876952214879858E10L,
        1.302170519734879977595901236693040544854E10L,
        2.922956950426397417800321486727032845006E9L,
        1.764317520783319397868923218385468729799E8L,
        1.573436014601118630105796794840834145120E7L,
        4.028077380105721388745632295157816229289E5L,
        1.644056806467289066852135096352853491530E4L,
        3.390868480059991640235675479463287886081E1L
    ];
    immutable real[10] TD1 = [
        -3.005357030696532927149885530689529032152E11L,
        -1.342602283126282827411658673839982164042E11L,
        -2.777153893355340961288511024443668743399E10L,
        -3.483826391033531996955620074072768276974E9L,
        -2.906321047071299585682722511260895227921E8L,
        -1.653347985722154162439387878512427542691E7L,
        -6.245520581562848778466500301865173123136E5L,
        -1.402124304177498828590239373389110545142E4L,
        -1.209368072473510674493129989468348633579E2L,
        1.0L
    ];
}
else
{
    /* erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
       1/128 <= 1/x < 1/8
       Peak relative error 1.9e-21  */
    immutable real[5] R = [ 0x1.b9f6d8b78e22459ep-6, 0x1.1b84686b0a4ea43ap-1,
       0x1.b8f6aebe96000c2ap+1, 0x1.cb1dbedac27c8ec2p+2, 0x1.cf885f8f572a4c14p+1
    ];

    immutable real[6] S = [
       0x1.87ae3cae5f65eb5ep-5, 0x1.01616f266f306d08p+0, 0x1.a4abe0411eed6c22p+2,
       0x1.eac9ce3da600abaap+3, 0x1.5752a9ac2faebbccp+3, 1.0
    ];

    /* erf(x)  = x P(x^2)/Q(x^2)
       0 <= x <= 1
       Peak relative error 7.6e-23  */
    immutable real[7] T = [ 0x1.0da01654d757888cp+20, 0x1.2eb7497bc8b4f4acp+17,
       0x1.79078c19530f72a8p+15, 0x1.4eaf2126c0b2c23p+11, 0x1.1f2ea81c9d272a2ep+8,
       0x1.59ca6e2d866e625p+2, 0x1.c188e0b67435faf4p-4
    ];

    immutable real[7] U = [ 0x1.dde6025c395ae34ep+19, 0x1.c4bc8b6235df35aap+18,
       0x1.8465900e88b6903ap+16, 0x1.855877093959ffdp+13, 0x1.e5c44395625ee358p+9,
       0x1.6a0fed103f1c68a6p+5, 1.0
    ];
}
}

/**
 *  Complementary error function
 *
 * erfc(x) = 1 - erf(x), and has high relative accuracy for
 * values of x far from zero. (For values near zero, use erf(x)).
 *
 *  1 - erf(x) =  2/ $(SQRT)(&pi;)
 *     $(INTEGRAL x, $(INFINITY)) exp( - $(POWER t, 2)) dt
 *
 *
 * For small x, erfc(x) = 1 - erf(x); otherwise rational
 * approximations are computed.
 *
 * A special function expx2(x) is used to suppress error amplification
 * in computing exp(-x^2).
 */
real erfc(real a)
{
    if (a == real.infinity)
        return 0.0;
    if (a == -real.infinity)
        return 2.0;

    immutable x = (a < 0.0) ? -a : a;

    if (x < (isIEEEQuadruple ? 0.25 : 1.0))
        return 1.0 - erf(a);

    static if (isIEEEQuadruple)
    {
        if (x < 1.25)
        {
            real y;
            final switch (cast(int)(8.0 * x))
            {
                case 2:
                    const z = x - 0.25;
                    y = C13b + z * rationalPoly(z, RNr13, RDr13);
                    y += C13a;
                    break;
                case 3:
                    const z = x - 0.375;
                    y = C14b + z * rationalPoly(z, RNr14, RDr14);
                    y += C14a;
                    break;
                case 4:
                    const z = x - 0.5;
                    y = C15b + z * rationalPoly(z, RNr15, RDr15);
                    y += C15a;
                    break;
                case 5:
                    const z = x - 0.625;
                    y = C16b + z * rationalPoly(z, RNr16, RDr16);
                    y += C16a;
                    break;
                case 6:
                    const z = x - 0.75;
                    y = C17b + z * rationalPoly(z, RNr17, RDr17);
                    y += C17a;
                    break;
                case 7:
                    const z = x - 0.875;
                    y = C18b + z * rationalPoly(z, RNr18, RDr18);
                    y += C18a;
                    break;
                case 8:
                    const z = x - 1.0;
                    y = C19b + z * rationalPoly(z, RNr19, RDr19);
                    y += C19a;
                    break;
                case 9:
                    const z = x - 1.125;
                    y = C20b + z * rationalPoly(z, RNr20, RDr20);
                    y += C20a;
                    break;
            }
            if (a < 0.0)
                y = 2.0 - y;
            return y;
        }
    }

    if (-a * a < -MAXLOG)
    {
        // underflow
        if (a < 0.0) return 2.0;
        else return 0.0;
    }

    real y;
    immutable z = expx2(a, -1);

    static if (isIEEEQuadruple)
    {
        y = z * erfce(x);
    }
    else
    {
        y = 1.0 / x;
        if (x < 8.0)
            y = z * rationalPoly(y, P, Q);
        else
            y = z * y * rationalPoly(y * y, R, S);
    }

    if (a < 0.0)
        y = 2.0 - y;

    if (y == 0.0)
    {
        // underflow
        if (a < 0.0) return 2.0;
        else return 0.0;
    }

    return y;
}


private {
/* Exponentially scaled erfc function
   exp(x^2) erfc(x)
   valid for x > 1.
   Use with normalDistribution and expx2.  */
static if (isIEEEQuadruple)
{
    real erfce(real x)
    {
        immutable z = 1.0L / (x * x);

        real p;
        switch (cast(int)(8.0 / x))
        {
            default:
            case 0:
                p = rationalPoly(z, RNr1, RDr1);
                break;
            case 1:
                p = rationalPoly(z, RNr2, RDr2);
                break;
            case 2:
                p = rationalPoly(z, RNr3, RDr3);
                break;
            case 3:
                p = rationalPoly(z, RNr4, RDr4);
                break;
            case 4:
                p = rationalPoly(z, RNr5, RDr5);
                break;
            case 5:
                p = rationalPoly(z, RNr6, RDr6);
                break;
            case 6:
                p = rationalPoly(z, RNr7, RDr7);
                break;
            case 7:
                p = rationalPoly(z, RNr8, RDr8);
                break;
        }
        return p / x;
    }
}
else
{
    real erfce(real x)
    {
        real y = 1.0/x;

        if (x < 8.0)
        {
            return rationalPoly(y, P, Q);
        }
        else
        {
            return y * rationalPoly(y * y, R, S);
        }
    }
}
}

/**
 *  Error function
 *
 * The integral is
 *
 *  erf(x) =  2/ $(SQRT)(&pi;)
 *     $(INTEGRAL 0, x) exp( - $(POWER t, 2)) dt
 *
 * The magnitude of x is limited to about 106.56 for IEEE 80-bit
 * arithmetic; 1 or -1 is returned outside this range.
 *
 * For 0 <= |x| < 1, a rational polynomials are used; otherwise
 * erf(x) = 1 - erfc(x).
 *
 * ACCURACY:
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0,1         50000       2.0e-19     5.7e-20
 */
real erf(real x)
{
    if (x == 0.0)
        return x; // deal with negative zero
    if (x == -real.infinity)
        return -1.0;
    if (x == real.infinity)
        return 1.0;
    immutable ax = abs(x);
    if (ax > 1.0L)
        return 1.0L - erfc(x);

    static if (isIEEEQuadruple)
    {
        immutable z = x * x;

        real y;
        if (ax < 0.875)
        {
            y = ax + ax * rationalPoly(x * x, TN1, TD1);
        }
        else
        {
            y = erfConst + rationalPoly(ax - 1.0L, TN2, TD2);
        }

        if (x < 0)
            y = -y;
        return y;
    }
    else
    {
        real z = x * x;
        return x * rationalPoly(x * x, T, U);
    }
}

@safe unittest
{
    // High resolution test points.
    enum real erfc0_250 = 0.723663330078125 + 1.0279753638067014931732235184287934646022E-5;
    enum real erfc0_375 = 0.5958709716796875 + 1.2118885490201676174914080878232469565953E-5;
    enum real erfc0_500 = 0.4794921875 + 7.9346869534623172533461080354712635484242E-6;
    enum real erfc0_625 = 0.3767547607421875 + 4.3570693945275513594941232097252997287766E-6;
    enum real erfc0_750 = 0.2888336181640625 + 1.0748182422368401062165408589222625794046E-5;
    enum real erfc0_875 = 0.215911865234375 + 1.3073705765341685464282101150637224028267E-5;
    enum real erfc1_000 = 0.15728759765625 + 1.1609394035130658779364917390740703933002E-5;
    enum real erfc1_125 = 0.111602783203125 + 8.9850951672359304215530728365232161564636E-6;

    enum real erf0_875  = (1-0.215911865234375) - 1.3073705765341685464282101150637224028267E-5;

    static bool isNaNWithPayload(real x, ulong payload) @safe pure nothrow @nogc
    {
      return isNaN(x) && getNaNPayload(x) == payload;
    }

    assert(feqrel(erfc(0.250L), erfc0_250 )>=real.mant_dig-1);
    assert(feqrel(erfc(0.375L), erfc0_375 )>=real.mant_dig-0);
    assert(feqrel(erfc(0.500L), erfc0_500 )>=real.mant_dig-2);
    assert(feqrel(erfc(0.625L), erfc0_625 )>=real.mant_dig-1);
    assert(feqrel(erfc(0.750L), erfc0_750 )>=real.mant_dig-1);
    assert(feqrel(erfc(0.875L), erfc0_875 )>=real.mant_dig-4);
    assert(feqrel(erfc(1.000L), erfc1_000 )>=real.mant_dig-2);
    assert(feqrel(erfc(1.125L), erfc1_125 )>=real.mant_dig-2);
    assert(feqrel(erf(0.875L), erf0_875 )>=real.mant_dig-1);
    // The DMC implementation of erfc() fails this next test (just).
    // Test point from Mathematica 11.0.
    assert(feqrel(erfc(4.1L), 6.70002765408489837272673380763418472e-9L) >= real.mant_dig-5);

    assert(isIdentical(erf(0.0),0.0));
    assert(isIdentical(erf(-0.0),-0.0));
    assert(erf(real.infinity) == 1.0);
    assert(erf(-real.infinity) == -1.0);
    assert(isNaNWithPayload(erf(NaN(0xDEF)), 0xDEF));
    assert(isNaNWithPayload(erfc(NaN(0xDEF)), 0xDEF));
    assert(isIdentical(erfc(real.infinity),0.0));
    assert(erfc(-real.infinity) == 2.0);
    assert(erfc(0) == 1.0);
}

/*
 *  Exponential of squared argument
 *
 * Computes y = exp(x*x) while suppressing error amplification
 * that would ordinarily arise from the inexactness of the
 * exponential argument x*x.
 *
 * If sign < 0, the result is inverted; i.e., y = exp(-x*x) .
 *
 * ACCURACY:
 *                      Relative error:
 * arithmetic      domain        # trials      peak         rms
 *   IEEE     -106.566, 106.566    10^5       1.6e-19     4.4e-20
 */

real expx2(real x, int sign)
{
    /*
    Cephes Math Library Release 2.9:  June, 2000
    Copyright 2000 by Stephen L. Moshier
    */
    const real M = 32_768.0;
    const real MINV = 3.0517578125e-5L;

    x = abs(x);
    if (sign < 0)
        x = -x;

  /* Represent x as an exact multiple of M plus a residual.
     M is a power of 2 chosen so that exp(m * m) does not overflow
     or underflow and so that |x - m| is small.  */
    real m = MINV * floor(M * x + 0.5L);
    real f = x - m;

    /* x^2 = m^2 + 2mf + f^2 */
    real u = m * m;
    real u1 = 2 * m * f  +  f * f;

    if (sign < 0)
    {
        u = -u;
        u1 = -u1;
    }

    if ((u+u1) > MAXLOG)
        return real.infinity;

    /* u is exact, u1 is small.  */
    return exp(u) * exp(u1);
}


/*
Computes the normal distribution function.

The normal (or Gaussian, or bell-shaped) distribution is
defined as:

normalDist(x) = 1/$(SQRT) &pi; $(INTEGRAL -$(INFINITY), x) exp( - $(POWER t, 2)/2) dt
    = 0.5 + 0.5 * erf(x/sqrt(2))
    = 0.5 * erfc(- x/sqrt(2))

To maintain accuracy at high values of x, use
normalDistribution(x) = 1 - normalDistribution(-x).

Accuracy:
Within a few bits of machine resolution over the entire
range.

References:
$(LINK http://www.netlib.org/cephes/ldoubdoc.html),
G. Marsaglia, "Evaluating the Normal Distribution",
Journal of Statistical Software <b>11</b>, (July 2004).
*/
real normalDistributionImpl(real a)
{
    real x = a * SQRT1_2;
    real z = abs(x);

    if ( z < 1.0 )
        return 0.5L + 0.5L * erf(x);
    else
    {
        real y = 0.5L * erfce(z);
        /* Multiply by exp(-x^2 / 2)  */
        z = expx2(a, -1);
        y = y * sqrt(z);
        if ( x > 0.0L )
            y = 1.0L - y;
        return y;
    }
}

@safe unittest
{
assert(fabs(normalDistributionImpl(1L) - (0.841344746068543))< 0.0000000000000005);
assert(isIdentical(normalDistributionImpl(NaN(0x325)), NaN(0x325)));
}

/*
 * Inverse of Normal distribution function
 *
 * Returns the argument, x, for which the area under the
 * Normal probability density function (integrated from
 * minus infinity to x) is equal to p.
 *
 * For small arguments 0 < p < exp(-2), the program computes
 * z = sqrt( -2 log(p) );  then the approximation is
 * x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z) .
 * For larger arguments,  x/sqrt(2 pi) = w + w^3 R(w^2)/S(w^2)) ,
 * where w = p - 0.5.
 */
// TODO: isIEEEQuadruple (128 bit) real implementation; not available from CEPHES.
real normalDistributionInvImpl(real p)
in {
  assert(p >= 0.0L && p <= 1.0L, "Domain error");
}
body
{
static immutable real[8] P0 =
[ -0x1.758f4d969484bfdcp-7, 0x1.53cee17a59259dd2p-3,
   -0x1.ea01e4400a9427a2p-1,  0x1.61f7504a0105341ap+1, -0x1.09475a594d0399f6p+2,
    0x1.7c59e7a0df99e3e2p+1, -0x1.87a81da52edcdf14p-1,  0x1.1fb149fd3f83600cp-7
];

static immutable real[8] Q0 =
[ -0x1.64b92ae791e64bb2p-7, 0x1.7585c7d597298286p-3,
   -0x1.40011be4f7591ce6p+0, 0x1.1fc067d8430a425ep+2, -0x1.21008ffb1e7ccdf2p+3,
   0x1.3d1581cf9bc12fccp+3, -0x1.53723a89fd8f083cp+2, 1.0
];

static immutable real[10] P1 =
[ 0x1.20ceea49ea142f12p-13, 0x1.cbe8a7267aea80bp-7,
   0x1.79fea765aa787c48p-2, 0x1.d1f59faa1f4c4864p+1, 0x1.1c22e426a013bb96p+4,
   0x1.a8675a0c51ef3202p+5, 0x1.75782c4f83614164p+6, 0x1.7a2f3d90948f1666p+6,
   0x1.5cd116ee4c088c3ap+5, 0x1.1361e3eb6e3cc20ap+2
];

static immutable real[10] Q1 =
[ 0x1.3a4ce1406cea98fap-13, 0x1.f45332623335cda2p-7,
   0x1.98f28bbd4b98db1p-2, 0x1.ec3b24f9c698091cp+1, 0x1.1cc56ecda7cf58e4p+4,
   0x1.92c6f7376bf8c058p+5, 0x1.4154c25aa47519b4p+6, 0x1.1b321d3b927849eap+6,
   0x1.403a5f5a4ce7b202p+4, 1.0
];

static immutable real[8] P2 =
[ 0x1.8c124a850116a6d8p-21, 0x1.534abda3c2fb90bap-13,
   0x1.29a055ec93a4718cp-7, 0x1.6468e98aad6dd474p-3, 0x1.3dab2ef4c67a601cp+0,
   0x1.e1fb3a1e70c67464p+1, 0x1.b6cce8035ff57b02p+2, 0x1.9f4c9e749ff35f62p+1
];

static immutable real[8] Q2 =
[ 0x1.af03f4fc0655e006p-21, 0x1.713192048d11fb2p-13,
   0x1.4357e5bbf5fef536p-7, 0x1.7fdac8749985d43cp-3, 0x1.4a080c813a2d8e84p+0,
   0x1.c3a4b423cdb41bdap+1, 0x1.8160694e24b5557ap+2, 1.0
];

static immutable real[8] P3 =
[ -0x1.55da447ae3806168p-34, -0x1.145635641f8778a6p-24,
 -0x1.abf46d6b48040128p-17, -0x1.7da550945da790fcp-11, -0x1.aa0b2a31157775fap-8,
   0x1.b11d97522eed26bcp-3, 0x1.1106d22f9ae89238p+1, 0x1.029a358e1e630f64p+1
];

static immutable real[8] Q3 =
[ -0x1.74022dd5523e6f84p-34, -0x1.2cb60d61e29ee836p-24,
   -0x1.d19e6ec03a85e556p-17, -0x1.9ea2a7b4422f6502p-11, -0x1.c54b1e852f107162p-8,
   0x1.e05268dd3c07989ep-3, 0x1.239c6aff14afbf82p+1, 1.0
];

    if (p <= 0.0L || p >= 1.0L)
    {
        if (p == 0.0L)
            return -real.infinity;
        if ( p == 1.0L )
            return real.infinity;
        return real.nan; // domain error
    }
    int code = 1;
    real y = p;
    if ( y > (1.0L - EXP_2) )
    {
        y = 1.0L - y;
        code = 0;
    }

    real x, z, y2, x0, x1;

    if ( y > EXP_2 )
    {
        y = y - 0.5L;
        y2 = y * y;
        x = y + y * (y2 * rationalPoly( y2, P0, Q0));
        return x * SQRT2PI;
    }

    x = sqrt( -2.0L * log(y) );
    x0 = x - log(x)/x;
    z = 1.0L/x;
    if ( x < 8.0L )
    {
        x1 = z * rationalPoly( z, P1, Q1);
    }
    else if ( x < 32.0L )
    {
        x1 = z * rationalPoly( z, P2, Q2);
    }
    else
    {
        x1 = z * rationalPoly( z, P3, Q3);
    }
    x = x0 - x1;
    if ( code != 0 )
    {
        x = -x;
    }
    return x;
}


@safe unittest
{
    // TODO: Use verified test points.
    // The values below are from Excel 2003.
    assert(fabs(normalDistributionInvImpl(0.001) - (-3.09023230616779))< 0.00000000000005);
    assert(fabs(normalDistributionInvImpl(1e-50) - (-14.9333375347885))< 0.00000000000005);
    assert(feqrel(normalDistributionInvImpl(0.999), -normalDistributionInvImpl(0.001)) > real.mant_dig-6);

    // Excel 2003 gets all the following values wrong!
    assert(normalDistributionInvImpl(0.0) == -real.infinity);
    assert(normalDistributionInvImpl(1.0) == real.infinity);
    assert(normalDistributionInvImpl(0.5) == 0);
    // (Excel 2003 returns norminv(p) = -30 for all p < 1e-200).
    // The value tested here is the one the function returned in Jan 2006.
    real unknown1 = normalDistributionInvImpl(1e-250L);
    assert( fabs(unknown1 -(-33.79958617269L) ) < 0.00000005);
}