/**
* 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 = Γ
* INTEGRAL = ∫
* INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>)
* POWER = $1<sup>$2</sup>
* BIGSUM = $(BIG Σ <sup>$2</sup><sub>$(SMALL $1)</sub>)
* CHOOSE = $(BIG () <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG ))
* 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)(π)
* $(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)(π)
* $(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) π $(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);
}