1 /**
2 * Error Functions and Normal Distribution.
3 *
4 * License: $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0).
5 * Copyright: Based on the CEPHES math library, which is
6 * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
7 * Authors: Stephen L. Moshier, ported to D by Don Clugston and David Nadlinger
8 */
9 /**
10 * Macros:
11 * NAN = $(RED NAN)
12 * SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
13 * GAMMA = Γ
14 * INTEGRAL = ∫
15 * INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>)
16 * POWER = $1<sup>$2</sup>
17 * BIGSUM = $(BIG Σ <sup>$2</sup><sub>$(SMALL $1)</sub>)
18 * CHOOSE = $(BIG () <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG ))
19 * TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
20 * <caption>Special Values</caption>
21 * $0</table>
22 * SVH = $(TR $(TH $1) $(TH $2))
23 * SV = $(TR $(TD $1) $(TD $2))
24 */
25 module std.internal.math.errorfunction;
26 import std.math;
27
28 pure:
29 nothrow:
30 @safe:
31 @nogc:
32
33 private {
34 immutable real EXP_2 = 0.135335283236612691893999494972484403L; /* exp(-2) */
35 enum real SQRT2PI = 2.50662827463100050241576528481104525L; // sqrt(2pi)
36
37
38 enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max)
39 enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min*real.epsilon) = log(smallest denormal)
40 }
41
rationalPoly(T)42 T rationalPoly(T)(T x, const(T) [] numerator, const(T) [] denominator) pure nothrow
43 {
44 return poly(x, numerator)/poly(x, denominator);
45 }
46
47
48 private {
49
50 /* erfc(x) = exp(-x^2) P(1/x)/Q(1/x)
51 1/8 <= 1/x <= 1
52 Peak relative error 5.8e-21 */
53 immutable real[10] P = [ -0x1.30dfa809b3cc6676p-17, 0x1.38637cd0913c0288p+18,
54 0x1.2f015e047b4476bp+22, 0x1.24726f46aa9ab08p+25, 0x1.64b13c6395dc9c26p+27,
55 0x1.294c93046ad55b5p+29, 0x1.5962a82f92576dap+30, 0x1.11a709299faba04ap+31,
56 0x1.11028065b087be46p+31, 0x1.0d8ef40735b097ep+30
57 ];
58
59 immutable real[11] Q = [ 0x1.14d8e2a72dec49f4p+19, 0x1.0c880ff467626e1p+23,
60 0x1.04417ef060b58996p+26, 0x1.404e61ba86df4ebap+28, 0x1.0f81887bc82b873ap+30,
61 0x1.4552a5e39fb49322p+31, 0x1.11779a0ceb2a01cep+32, 0x1.3544dd691b5b1d5cp+32,
62 0x1.a91781f12251f02ep+31, 0x1.0d8ef3da605a1c86p+30, 1.0
63 ];
64
65 // For 128 bit quadruple-precision floats, we use a higher-precision implementation
66 // with more polynomial segments.
67 enum isIEEEQuadruple = floatTraits!real.realFormat == RealFormat.ieeeQuadruple;
68 static if (isIEEEQuadruple)
69 {
70 // erfc(x + 0.25) = erfc(0.25) + x R(x)
71 // 0 <= x < 0.125
72 // Peak relative error 1.4e-35
73 immutable real[9] RNr13 = [
74 -2.353707097641280550282633036456457014829E3L,
75 3.871159656228743599994116143079870279866E2L,
76 -3.888105134258266192210485617504098426679E2L,
77 -2.129998539120061668038806696199343094971E1L,
78 -8.125462263594034672468446317145384108734E1L,
79 8.151549093983505810118308635926270319660E0L,
80 -5.033362032729207310462422357772568553670E0L,
81 -4.253956621135136090295893547735851168471E-2L,
82 -8.098602878463854789780108161581050357814E-2L
83 ];
84 immutable real[9] RDr13 = [
85 2.220448796306693503549505450626652881752E3L,
86 1.899133258779578688791041599040951431383E2L,
87 1.061906712284961110196427571557149268454E3L,
88 7.497086072306967965180978101974566760042E1L,
89 2.146796115662672795876463568170441327274E2L,
90 1.120156008362573736664338015952284925592E1L,
91 2.211014952075052616409845051695042741074E1L,
92 6.469655675326150785692908453094054988938E-1L,
93 1.0
94 ];
95
96 // erfc(0.25) = C13a + C13b to extra precision.
97 immutable real C13a = 0.723663330078125L;
98 immutable real C13b = 1.0279753638067014931732235184287934646022E-5L;
99
100 // erfc(x + 0.375) = erfc(0.375) + x R(x)
101 // 0 <= x < 0.125
102 // Peak relative error 1.2e-35
103 immutable real[9] RNr14 = [
104 -2.446164016404426277577283038988918202456E3L,
105 6.718753324496563913392217011618096698140E2L,
106 -4.581631138049836157425391886957389240794E2L,
107 -2.382844088987092233033215402335026078208E1L,
108 -7.119237852400600507927038680970936336458E1L,
109 1.313609646108420136332418282286454287146E1L,
110 -6.188608702082264389155862490056401365834E0L,
111 -2.787116601106678287277373011101132659279E-2L,
112 -2.230395570574153963203348263549700967918E-2L
113 ];
114 immutable real[9] RDr14 = [
115 2.495187439241869732696223349840963702875E3L,
116 2.503549449872925580011284635695738412162E2L,
117 1.159033560988895481698051531263861842461E3L,
118 9.493751466542304491261487998684383688622E1L,
119 2.276214929562354328261422263078480321204E2L,
120 1.367697521219069280358984081407807931847E1L,
121 2.276988395995528495055594829206582732682E1L,
122 7.647745753648996559837591812375456641163E-1L,
123 1.0
124 ];
125
126 // erfc(0.375) = C14a + C14b to extra precision.
127 immutable real C14a = 0.5958709716796875L;
128 immutable real C14b = 1.2118885490201676174914080878232469565953E-5L;
129
130 // erfc(x + 0.5) = erfc(0.5) + x R(x)
131 // 0 <= x < 0.125
132 // Peak relative error 4.7e-36
133 immutable real[9] RNr15 = [
134 -2.624212418011181487924855581955853461925E3L,
135 8.473828904647825181073831556439301342756E2L,
136 -5.286207458628380765099405359607331669027E2L,
137 -3.895781234155315729088407259045269652318E1L,
138 -6.200857908065163618041240848728398496256E1L,
139 1.469324610346924001393137895116129204737E1L,
140 -6.961356525370658572800674953305625578903E0L,
141 5.145724386641163809595512876629030548495E-3L,
142 1.990253655948179713415957791776180406812E-2L
143 ];
144 immutable real[9] RDr15 = [
145 2.986190760847974943034021764693341524962E3L,
146 5.288262758961073066335410218650047725985E2L,
147 1.363649178071006978355113026427856008978E3L,
148 1.921707975649915894241864988942255320833E2L,
149 2.588651100651029023069013885900085533226E2L,
150 2.628752920321455606558942309396855629459E1L,
151 2.455649035885114308978333741080991380610E1L,
152 1.378826653595128464383127836412100939126E0L,
153 1.0
154 ];
155 // erfc(0.5) = C15a + C15b to extra precision.
156 immutable real C15a = 0.4794921875L;
157 immutable real C15b = 7.9346869534623172533461080354712635484242E-6L;
158
159 // erfc(x + 0.625) = erfc(0.625) + x R(x)
160 // 0 <= x < 0.125
161 // Peak relative error 5.1e-36
162 immutable real[9] RNr16 = [
163 -2.347887943200680563784690094002722906820E3L,
164 8.008590660692105004780722726421020136482E2L,
165 -5.257363310384119728760181252132311447963E2L,
166 -4.471737717857801230450290232600243795637E1L,
167 -4.849540386452573306708795324759300320304E1L,
168 1.140885264677134679275986782978655952843E1L,
169 -6.731591085460269447926746876983786152300E0L,
170 1.370831653033047440345050025876085121231E-1L,
171 2.022958279982138755020825717073966576670E-2L,
172 ];
173 immutable real[9] RDr16 = [
174 3.075166170024837215399323264868308087281E3L,
175 8.730468942160798031608053127270430036627E2L,
176 1.458472799166340479742581949088453244767E3L,
177 3.230423687568019709453130785873540386217E2L,
178 2.804009872719893612081109617983169474655E2L,
179 4.465334221323222943418085830026979293091E1L,
180 2.612723259683205928103787842214809134746E1L,
181 2.341526751185244109722204018543276124997E0L,
182 1.0
183 ];
184 // erfc(0.625) = C16a + C16b to extra precision.
185 immutable real C16a = 0.3767547607421875L;
186 immutable real C16b = 4.3570693945275513594941232097252997287766E-6L;
187
188 // erfc(x + 0.75) = erfc(0.75) + x R(x)
189 // 0 <= x < 0.125
190 // Peak relative error 1.7e-35
191 immutable real[9] RNr17 = [
192 -1.767068734220277728233364375724380366826E3L,
193 6.693746645665242832426891888805363898707E2L,
194 -4.746224241837275958126060307406616817753E2L,
195 -2.274160637728782675145666064841883803196E1L,
196 -3.541232266140939050094370552538987982637E1L,
197 6.988950514747052676394491563585179503865E0L,
198 -5.807687216836540830881352383529281215100E0L,
199 3.631915988567346438830283503729569443642E-1L,
200 -1.488945487149634820537348176770282391202E-2L
201 ];
202 immutable real[9] RDr17 = [
203 2.748457523498150741964464942246913394647E3L,
204 1.020213390713477686776037331757871252652E3L,
205 1.388857635935432621972601695296561952738E3L,
206 3.903363681143817750895999579637315491087E2L,
207 2.784568344378139499217928969529219886578E2L,
208 5.555800830216764702779238020065345401144E1L,
209 2.646215470959050279430447295801291168941E1L,
210 2.984905282103517497081766758550112011265E0L,
211 1.0
212 ];
213 // erfc(0.75) = C17a + C17b to extra precision.
214 immutable real C17a = 0.2888336181640625L;
215 immutable real C17b = 1.0748182422368401062165408589222625794046E-5L;
216
217
218 // erfc(x + 0.875) = erfc(0.875) + x R(x)
219 // 0 <= x < 0.125
220 // Peak relative error 2.2e-35
221 immutable real[9] RNr18 = [
222 -1.342044899087593397419622771847219619588E3L,
223 6.127221294229172997509252330961641850598E2L,
224 -4.519821356522291185621206350470820610727E2L,
225 1.223275177825128732497510264197915160235E1L,
226 -2.730789571382971355625020710543532867692E1L,
227 4.045181204921538886880171727755445395862E0L,
228 -4.925146477876592723401384464691452700539E0L,
229 5.933878036611279244654299924101068088582E-1L,
230 -5.557645435858916025452563379795159124753E-2L
231 ];
232 immutable real[9] RDr18 = [
233 2.557518000661700588758505116291983092951E3L,
234 1.070171433382888994954602511991940418588E3L,
235 1.344842834423493081054489613250688918709E3L,
236 4.161144478449381901208660598266288188426E2L,
237 2.763670252219855198052378138756906980422E2L,
238 5.998153487868943708236273854747564557632E1L,
239 2.657695108438628847733050476209037025318E1L,
240 3.252140524394421868923289114410336976512E0L,
241 1.0
242 ];
243
244 // erfc(0.875) = C18a + C18b to extra precision.
245 immutable real C18a = 0.215911865234375L;
246 immutable real C18b = 1.3073705765341685464282101150637224028267E-5L;
247
248 // erfc(x + 1.0) = erfc(1.0) + x R(x)
249 // 0 <= x < 0.125
250 // Peak relative error 1.6e-35
251 immutable real[9] RNr19 = [
252 -1.139180936454157193495882956565663294826E3L,
253 6.134903129086899737514712477207945973616E2L,
254 -4.628909024715329562325555164720732868263E2L,
255 4.165702387210732352564932347500364010833E1L,
256 -2.286979913515229747204101330405771801610E1L,
257 1.870695256449872743066783202326943667722E0L,
258 -4.177486601273105752879868187237000032364E0L,
259 7.533980372789646140112424811291782526263E-1L,
260 -8.629945436917752003058064731308767664446E-2L
261 ];
262 immutable real[9] RDr19 = [
263 2.744303447981132701432716278363418643778E3L,
264 1.266396359526187065222528050591302171471E3L,
265 1.466739461422073351497972255511919814273E3L,
266 4.868710570759693955597496520298058147162E2L,
267 2.993694301559756046478189634131722579643E2L,
268 6.868976819510254139741559102693828237440E1L,
269 2.801505816247677193480190483913753613630E1L,
270 3.604439909194350263552750347742663954481E0L,
271 1.0
272 ];
273
274 // erfc(1.0) = C19a + C19b to extra precision.
275 immutable real C19a = 0.15728759765625L;
276 immutable real C19b = 1.1609394035130658779364917390740703933002E-5L;
277
278 // erfc(x + 1.125) = erfc(1.125) + x R(x)
279 // 0 <= x < 0.125
280 // Peak relative error 3.6e-36
281 immutable real[9] RNr20 = [
282 -9.652706916457973956366721379612508047640E2L,
283 5.577066396050932776683469951773643880634E2L,
284 -4.406335508848496713572223098693575485978E2L,
285 5.202893466490242733570232680736966655434E1L,
286 -1.931311847665757913322495948705563937159E1L,
287 -9.364318268748287664267341457164918090611E-2L,
288 -3.306390351286352764891355375882586201069E0L,
289 7.573806045289044647727613003096916516475E-1L,
290 -9.611744011489092894027478899545635991213E-2L
291 ];
292 immutable real[9] RDr20 = [
293 3.032829629520142564106649167182428189014E3L,
294 1.659648470721967719961167083684972196891E3L,
295 1.703545128657284619402511356932569292535E3L,
296 6.393465677731598872500200253155257708763E2L,
297 3.489131397281030947405287112726059221934E2L,
298 8.848641738570783406484348434387611713070E1L,
299 3.132269062552392974833215844236160958502E1L,
300 4.430131663290563523933419966185230513168E0L,
301 1.0
302 ];
303
304 // erfc(1.125) = C20a + C20b to extra precision.
305 immutable real C20a = 0.111602783203125L;
306 immutable real C20b = 8.9850951672359304215530728365232161564636E-6L;
307
308 // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
309 // 7/8 <= 1/x < 1
310 // Peak relative error 1.4e-35
311 immutable real[10] RNr8 = [
312 3.587451489255356250759834295199296936784E1L,
313 5.406249749087340431871378009874875889602E2L,
314 2.931301290625250886238822286506381194157E3L,
315 7.359254185241795584113047248898753470923E3L,
316 9.201031849810636104112101947312492532314E3L,
317 5.749697096193191467751650366613289284777E3L,
318 1.710415234419860825710780802678697889231E3L,
319 2.150753982543378580859546706243022719599E2L,
320 8.740953582272147335100537849981160931197E0L,
321 4.876422978828717219629814794707963640913E-2L
322 ];
323 immutable real[10] RDr8 = [
324 6.358593134096908350929496535931630140282E1L,
325 9.900253816552450073757174323424051765523E2L,
326 5.642928777856801020545245437089490805186E3L,
327 1.524195375199570868195152698617273739609E4L,
328 2.113829644500006749947332935305800887345E4L,
329 1.526438562626465706267943737310282977138E4L,
330 5.561370922149241457131421914140039411782E3L,
331 9.394035530179705051609070428036834496942E2L,
332 6.147019596150394577984175188032707343615E1L,
333 1.0L
334 ];
335
336 // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
337 // 3/4 <= 1/x < 7/8
338 // Peak relative error 1.7e-36
339 immutable real[10] RNr7 = [
340 1.293515364043117601705535485785956592493E2L,
341 2.474534371269189867053251150063459055230E3L,
342 1.756537563959875738809491329250457486510E4L,
343 5.977479535376639763773153344676726091607E4L,
344 1.054313816139671870123172936972055385049E5L,
345 9.754699773487726957401038094714603033904E4L,
346 4.579131242577671038339922925213209214880E4L,
347 1.000710322164510887997115157797717324370E4L,
348 8.496863250712471449526805271633794700452E2L,
349 1.797349831386892396933210199236530557333E1L
350 ];
351 immutable real[11] RDr7 = [
352 2.292696320307033494820399866075534515002E2L,
353 4.500632608295626968062258401895610053116E3L,
354 3.321218723485498111535866988511716659339E4L,
355 1.196084512221845156596781258490840961462E5L,
356 2.287033883912529843927983406878910939930E5L,
357 2.370223495794642027268482075021298394425E5L,
358 1.305173734022437154610938308944995159199E5L,
359 3.589386258485887630236490009835928559621E4L,
360 4.339996864041074149726360516336773136101E3L,
361 1.753135522665469574605384979152863899099E2L,
362 1.0L
363 ];
364
365 // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
366 // 5/8 <= 1/x < 3/4
367 // Peak relative error 1.6e-35
368 immutable real[10] RNr6 = [
369 1.423313561367394923305025174137639124533E1L,
370 3.244462503273003837514629113846075327206E2L,
371 2.784937282403293364911673341412846781934E3L,
372 1.163060685810874867196849890286455473382E4L,
373 2.554141394931962276102205517358731053756E4L,
374 2.982733782500729530503336931258698708782E4L,
375 1.789683564523810605328169719436374742840E4L,
376 5.056032142227470121262177112822018882754E3L,
377 5.605349942234782054561269306895707034586E2L,
378 1.561652599080729507274832243665726064881E1L
379 ];
380 immutable real[11] RDr6 = [
381 2.522757606611479946069351519410222913326E1L,
382 5.876797910931896554014229647006604017806E2L,
383 5.211092128250480712011248211246144751074E3L,
384 2.282679910855404599271496827409168580797E4L,
385 5.371245819205596609986320599133109262447E4L,
386 6.926186102106400355114925675028888924445E4L,
387 4.794366033363621432575096487724913414473E4L,
388 1.673190682734065914573814938835674963896E4L,
389 2.589544846151313120096957014256536236242E3L,
390 1.349438432583208276883323156200117027433E2L,
391 1.0L
392 ];
393
394 // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
395 // 1/2 <= 1/x < 5/8
396 // Peak relative error 4.3e-36
397 immutable real[11] RNr5 = [
398 6.743447478279267339131211137241149796763E-2L,
399 2.031339015932962998168472743355874796350E0L,
400 2.369234815713196480221800940201367484379E1L,
401 1.387819614016107433603101545594790875922E2L,
402 4.435600256936515839937720907171966121786E2L,
403 7.881577949936817507981170892417739733047E2L,
404 7.615749099291545976179905281851765734680E2L,
405 3.752484528663442467089606663006771157777E2L,
406 8.279644286027145214308303292537009564726E1L,
407 6.201462983413738162709722770960040042647E0L,
408 6.649631608720062333043506249503378282697E-2L
409 ];
410 immutable real[11] RDr5 = [
411 1.195244945161889822018178270706903972343E-1L,
412 3.660216908153253021384862427197665991311E0L,
413 4.373405883243078019655721779021995159854E1L,
414 2.653305963056235008916733402786877121865E2L,
415 8.921329790491152046318422124415895506335E2L,
416 1.705552231555600759729260640562363304312E3L,
417 1.832711109606893446763174603477244625325E3L,
418 1.056823953275835649973998168744261083316E3L,
419 2.975561981792909722126456997074344895584E2L,
420 3.393149095158232521894537008472203487436E1L,
421 1.0L
422 ];
423
424 // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
425 // 3/8 <= 1/x < 1/2
426 // Peak relative error 1.8e-36
427 immutable real[11] RNr4 = [
428 3.558685919236420073872459554885612994007E-2L,
429 1.460223642496950651561817195253277924528E0L,
430 2.379856746555189546876720308841066577268E1L,
431 2.005205521246554860334064698817220160117E2L,
432 9.533374928664989955629120027419490517596E2L,
433 2.623024576994438336130421711314560425373E3L,
434 4.126446434603735586340585027628851620886E3L,
435 3.540675861596687801829655387867654030013E3L,
436 1.506037084891064572653273761987617394259E3L,
437 2.630715699182706745867272452228891752353E2L,
438 1.202476629652900619635409242749750364878E1L
439 ];
440 immutable real[12] RDr4 = [
441 6.307606561714590590399683184410336583739E-2L,
442 2.619717051134271249293056836082721776665E0L,
443 4.344441402681725017630451522968410844608E1L,
444 3.752891116408399440953195184301023399176E2L,
445 1.849305988804654653921972804388006355502E3L,
446 5.358505261991675891835885654499883449403E3L,
447 9.091890995405251314631428721090705475825E3L,
448 8.731418313949291797856351745278287516416E3L,
449 4.420211285043270337492325400764271868740E3L,
450 1.031487363021856106882306509107923200832E3L,
451 8.387036084846046121805145056040429461783E1L,
452 1.0L
453 ];
454
455 // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
456 // 1/4 <= 1/x < 3/8
457 // Peak relative error 8.1e-37
458 immutable real[12] RNr3 = [
459 4.584481542956275354582319313040418316755E-5L,
460 2.674923158288848442110883948437930349128E-3L,
461 6.344232532055212248017211243740466847311E-2L,
462 7.985145965992002744933550450451513513963E-1L,
463 5.845078061888281665064746347663123946270E0L,
464 2.566625318816866587482759497608029522596E1L,
465 6.736225182343446605268837827950856640948E1L,
466 1.021796460139598089409347761712730512053E2L,
467 8.344336615515430530929955615400706619764E1L,
468 3.207749011528249356283897356277376306967E1L,
469 4.386185123863412086856423971695142026036E0L,
470 8.971576448581208351826868348023528863856E-2L
471 ];
472 immutable real[12] RDr3 = [
473 8.125781965218112303281657065320409661370E-5L,
474 4.781806762611504685247817818428945295520E-3L,
475 1.147785538413798317790357996845767614561E-1L,
476 1.469285552007088106614218996464752307606E0L,
477 1.101712261349880339221039938999124077650E1L,
478 5.008507527095093413713171655268276861071E1L,
479 1.383058691613468970486425146336829447433E2L,
480 2.264114250278912520501010108736339599752E2L,
481 2.081377197698598680576330179979996940039E2L,
482 9.724438129690013609440151781601781137944E1L,
483 1.907905050771832372089975877589291760121E1L,
484 1.0L
485 ];
486
487 // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
488 // 1/8 <= 1/x < 1/4
489 // Peak relative error 1.5e-36
490 immutable real[12] RNr2 = [
491 6.928615158005256885698045840589513728399E-7L,
492 5.616245938942075826026382337922413007879E-5L,
493 1.871624980715261794832358438894219696113E-3L,
494 3.349922063795792371642023765253747563009E-2L,
495 3.531865233974349943956345502463135695834E-1L,
496 2.264714157625072773976468825160906342360E0L,
497 8.810720294489253776747319730638214883026E0L,
498 2.014056685571655833019183248931442888437E1L,
499 2.524586947657190747039554310814128743320E1L,
500 1.520656940937208886246188940244581671609E1L,
501 3.334145500790963675035841482334493680498E0L,
502 1.122108380007109245896534245151140632457E-1L
503 ];
504 immutable real[12] RDr2 = [
505 1.228065061824874795984937092427781089256E-6L,
506 1.001593999520159167559129042893802235969E-4L,
507 3.366527555699367241421450749821030974446E-3L,
508 6.098626947195865254152265585991861150369E-2L,
509 6.541547922508613985813189387198804660235E-1L,
510 4.301130233305371976727117480925676583204E0L,
511 1.737155892350891711527711121692994762909E1L,
512 4.206892112110558214680649401236873828801E1L,
513 5.787487996025016843403524261574779631219E1L,
514 4.094047148590822715163181507813774861621E1L,
515 1.230603930056944875836549716515643997094E1L,
516 1.0L
517 ];
518
519 // erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
520 // 1/128 <= 1/x < 1/8
521 // Peak relative error 2.2e-36
522 immutable real[10] RNr1 = [
523 1.293111801138199795250035229383033322539E-6L,
524 9.785224720880980456171438759161161816706E-5L,
525 2.932474396233212166056331430621176065943E-3L,
526 4.496429595719847083917435337780697436921E-2L,
527 3.805989855927720844877478869846718877846E-1L,
528 1.789745532222426292126781724570152590071E0L,
529 4.465737379634389318903237306594171764628E0L,
530 5.268535822258082278401240171488850433767E0L,
531 2.258357766807433839494276681092713991651E0L,
532 1.504459334078750002966538036652860809497E-1L
533 ];
534 immutable real[10] RDr1 = [
535 2.291980991578770070179177302906728371406E-6L,
536 1.745845828808028552029674694534934620384E-4L,
537 5.283248841982102317072923869576785278019E-3L,
538 8.221212297078141470944454807434634848018E-2L,
539 7.120500674861902950423510939060230945621E-1L,
540 3.475435367560809622183983439133664598155E0L,
541 9.243253391989233533874386043611304387113E0L,
542 1.227894792475280941511758877318903197188E1L,
543 6.789361410398841316638617624392719077724E0L,
544 1.0L
545 ];
546
547 // erf(z+1) = erfConst + P(z)/Q(z)
548 // -.125 <= z <= 0
549 // Peak relative error 7.3e-36
550 immutable real erfConst = 0.845062911510467529296875L;
551 immutable real[9] TN2 = [
552 -4.088889697077485301010486931817357000235E1L,
553 7.157046430681808553842307502826960051036E3L,
554 -2.191561912574409865550015485451373731780E3L,
555 2.180174916555316874988981177654057337219E3L,
556 2.848578658049670668231333682379720943455E2L,
557 1.630362490952512836762810462174798925274E2L,
558 6.317712353961866974143739396865293596895E0L,
559 2.450441034183492434655586496522857578066E1L,
560 5.127662277706787664956025545897050896203E-1L
561 ];
562 immutable real[10] TD2 = [
563 1.731026445926834008273768924015161048885E4L,
564 1.209682239007990370796112604286048173750E4L,
565 1.160950290217993641320602282462976163857E4L,
566 5.394294645127126577825507169061355698157E3L,
567 2.791239340533632669442158497532521776093E3L,
568 8.989365571337319032943005387378993827684E2L,
569 2.974016493766349409725385710897298069677E2L,
570 6.148192754590376378740261072533527271947E1L,
571 1.178502892490738445655468927408440847480E1L,
572 1.0L
573 ];
574
575 // erf(x) = x + x P(x^2)/Q(x^2)
576 // 0 <= x <= 7/8
577 // Peak relative error 1.8e-35
578 immutable real[9] TN1 = [
579 -3.858252324254637124543172907442106422373E10L,
580 9.580319248590464682316366876952214879858E10L,
581 1.302170519734879977595901236693040544854E10L,
582 2.922956950426397417800321486727032845006E9L,
583 1.764317520783319397868923218385468729799E8L,
584 1.573436014601118630105796794840834145120E7L,
585 4.028077380105721388745632295157816229289E5L,
586 1.644056806467289066852135096352853491530E4L,
587 3.390868480059991640235675479463287886081E1L
588 ];
589 immutable real[10] TD1 = [
590 -3.005357030696532927149885530689529032152E11L,
591 -1.342602283126282827411658673839982164042E11L,
592 -2.777153893355340961288511024443668743399E10L,
593 -3.483826391033531996955620074072768276974E9L,
594 -2.906321047071299585682722511260895227921E8L,
595 -1.653347985722154162439387878512427542691E7L,
596 -6.245520581562848778466500301865173123136E5L,
597 -1.402124304177498828590239373389110545142E4L,
598 -1.209368072473510674493129989468348633579E2L,
599 1.0L
600 ];
601 }
602 else
603 {
604 /* erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
605 1/128 <= 1/x < 1/8
606 Peak relative error 1.9e-21 */
607 immutable real[5] R = [ 0x1.b9f6d8b78e22459ep-6, 0x1.1b84686b0a4ea43ap-1,
608 0x1.b8f6aebe96000c2ap+1, 0x1.cb1dbedac27c8ec2p+2, 0x1.cf885f8f572a4c14p+1
609 ];
610
611 immutable real[6] S = [
612 0x1.87ae3cae5f65eb5ep-5, 0x1.01616f266f306d08p+0, 0x1.a4abe0411eed6c22p+2,
613 0x1.eac9ce3da600abaap+3, 0x1.5752a9ac2faebbccp+3, 1.0
614 ];
615
616 /* erf(x) = x P(x^2)/Q(x^2)
617 0 <= x <= 1
618 Peak relative error 7.6e-23 */
619 immutable real[7] T = [ 0x1.0da01654d757888cp+20, 0x1.2eb7497bc8b4f4acp+17,
620 0x1.79078c19530f72a8p+15, 0x1.4eaf2126c0b2c23p+11, 0x1.1f2ea81c9d272a2ep+8,
621 0x1.59ca6e2d866e625p+2, 0x1.c188e0b67435faf4p-4
622 ];
623
624 immutable real[7] U = [ 0x1.dde6025c395ae34ep+19, 0x1.c4bc8b6235df35aap+18,
625 0x1.8465900e88b6903ap+16, 0x1.855877093959ffdp+13, 0x1.e5c44395625ee358p+9,
626 0x1.6a0fed103f1c68a6p+5, 1.0
627 ];
628 }
629 }
630
631 /**
632 * Complementary error function
633 *
634 * erfc(x) = 1 - erf(x), and has high relative accuracy for
635 * values of x far from zero. (For values near zero, use erf(x)).
636 *
637 * 1 - erf(x) = 2/ $(SQRT)(π)
638 * $(INTEGRAL x, $(INFINITY)) exp( - $(POWER t, 2)) dt
639 *
640 *
641 * For small x, erfc(x) = 1 - erf(x); otherwise rational
642 * approximations are computed.
643 *
644 * A special function expx2(x) is used to suppress error amplification
645 * in computing exp(-x^2).
646 */
erfc(real a)647 real erfc(real a)
648 {
649 if (a == real.infinity)
650 return 0.0;
651 if (a == -real.infinity)
652 return 2.0;
653
654 immutable x = (a < 0.0) ? -a : a;
655
656 if (x < (isIEEEQuadruple ? 0.25 : 1.0))
657 return 1.0 - erf(a);
658
659 static if (isIEEEQuadruple)
660 {
661 if (x < 1.25)
662 {
663 real y;
664 final switch (cast(int)(8.0 * x))
665 {
666 case 2:
667 const z = x - 0.25;
668 y = C13b + z * rationalPoly(z, RNr13, RDr13);
669 y += C13a;
670 break;
671 case 3:
672 const z = x - 0.375;
673 y = C14b + z * rationalPoly(z, RNr14, RDr14);
674 y += C14a;
675 break;
676 case 4:
677 const z = x - 0.5;
678 y = C15b + z * rationalPoly(z, RNr15, RDr15);
679 y += C15a;
680 break;
681 case 5:
682 const z = x - 0.625;
683 y = C16b + z * rationalPoly(z, RNr16, RDr16);
684 y += C16a;
685 break;
686 case 6:
687 const z = x - 0.75;
688 y = C17b + z * rationalPoly(z, RNr17, RDr17);
689 y += C17a;
690 break;
691 case 7:
692 const z = x - 0.875;
693 y = C18b + z * rationalPoly(z, RNr18, RDr18);
694 y += C18a;
695 break;
696 case 8:
697 const z = x - 1.0;
698 y = C19b + z * rationalPoly(z, RNr19, RDr19);
699 y += C19a;
700 break;
701 case 9:
702 const z = x - 1.125;
703 y = C20b + z * rationalPoly(z, RNr20, RDr20);
704 y += C20a;
705 break;
706 }
707 if (a < 0.0)
708 y = 2.0 - y;
709 return y;
710 }
711 }
712
713 if (-a * a < -MAXLOG)
714 {
715 // underflow
716 if (a < 0.0) return 2.0;
717 else return 0.0;
718 }
719
720 real y;
721 immutable z = expx2(a, -1);
722
723 static if (isIEEEQuadruple)
724 {
725 y = z * erfce(x);
726 }
727 else
728 {
729 y = 1.0 / x;
730 if (x < 8.0)
731 y = z * rationalPoly(y, P, Q);
732 else
733 y = z * y * rationalPoly(y * y, R, S);
734 }
735
736 if (a < 0.0)
737 y = 2.0 - y;
738
739 if (y == 0.0)
740 {
741 // underflow
742 if (a < 0.0) return 2.0;
743 else return 0.0;
744 }
745
746 return y;
747 }
748
749
750 private {
751 /* Exponentially scaled erfc function
752 exp(x^2) erfc(x)
753 valid for x > 1.
754 Use with normalDistribution and expx2. */
755 static if (isIEEEQuadruple)
756 {
erfce(real x)757 real erfce(real x)
758 {
759 immutable z = 1.0L / (x * x);
760
761 real p;
762 switch (cast(int)(8.0 / x))
763 {
764 default:
765 case 0:
766 p = rationalPoly(z, RNr1, RDr1);
767 break;
768 case 1:
769 p = rationalPoly(z, RNr2, RDr2);
770 break;
771 case 2:
772 p = rationalPoly(z, RNr3, RDr3);
773 break;
774 case 3:
775 p = rationalPoly(z, RNr4, RDr4);
776 break;
777 case 4:
778 p = rationalPoly(z, RNr5, RDr5);
779 break;
780 case 5:
781 p = rationalPoly(z, RNr6, RDr6);
782 break;
783 case 6:
784 p = rationalPoly(z, RNr7, RDr7);
785 break;
786 case 7:
787 p = rationalPoly(z, RNr8, RDr8);
788 break;
789 }
790 return p / x;
791 }
792 }
793 else
794 {
erfce(real x)795 real erfce(real x)
796 {
797 real y = 1.0/x;
798
799 if (x < 8.0)
800 {
801 return rationalPoly(y, P, Q);
802 }
803 else
804 {
805 return y * rationalPoly(y * y, R, S);
806 }
807 }
808 }
809 }
810
811 /**
812 * Error function
813 *
814 * The integral is
815 *
816 * erf(x) = 2/ $(SQRT)(π)
817 * $(INTEGRAL 0, x) exp( - $(POWER t, 2)) dt
818 *
819 * The magnitude of x is limited to about 106.56 for IEEE 80-bit
820 * arithmetic; 1 or -1 is returned outside this range.
821 *
822 * For 0 <= |x| < 1, a rational polynomials are used; otherwise
823 * erf(x) = 1 - erfc(x).
824 *
825 * ACCURACY:
826 * Relative error:
827 * arithmetic domain # trials peak rms
828 * IEEE 0,1 50000 2.0e-19 5.7e-20
829 */
erf(real x)830 real erf(real x)
831 {
832 if (x == 0.0)
833 return x; // deal with negative zero
834 if (x == -real.infinity)
835 return -1.0;
836 if (x == real.infinity)
837 return 1.0;
838 immutable ax = abs(x);
839 if (ax > 1.0L)
840 return 1.0L - erfc(x);
841
842 static if (isIEEEQuadruple)
843 {
844 immutable z = x * x;
845
846 real y;
847 if (ax < 0.875)
848 {
849 y = ax + ax * rationalPoly(x * x, TN1, TD1);
850 }
851 else
852 {
853 y = erfConst + rationalPoly(ax - 1.0L, TN2, TD2);
854 }
855
856 if (x < 0)
857 y = -y;
858 return y;
859 }
860 else
861 {
862 real z = x * x;
863 return x * rationalPoly(x * x, T, U);
864 }
865 }
866
867 @safe unittest
868 {
869 // High resolution test points.
870 enum real erfc0_250 = 0.723663330078125 + 1.0279753638067014931732235184287934646022E-5;
871 enum real erfc0_375 = 0.5958709716796875 + 1.2118885490201676174914080878232469565953E-5;
872 enum real erfc0_500 = 0.4794921875 + 7.9346869534623172533461080354712635484242E-6;
873 enum real erfc0_625 = 0.3767547607421875 + 4.3570693945275513594941232097252997287766E-6;
874 enum real erfc0_750 = 0.2888336181640625 + 1.0748182422368401062165408589222625794046E-5;
875 enum real erfc0_875 = 0.215911865234375 + 1.3073705765341685464282101150637224028267E-5;
876 enum real erfc1_000 = 0.15728759765625 + 1.1609394035130658779364917390740703933002E-5;
877 enum real erfc1_125 = 0.111602783203125 + 8.9850951672359304215530728365232161564636E-6;
878
879 enum real erf0_875 = (1-0.215911865234375) - 1.3073705765341685464282101150637224028267E-5;
880
isNaNWithPayload(real x,ulong payload)881 static bool isNaNWithPayload(real x, ulong payload) @safe pure nothrow @nogc
882 {
883 return isNaN(x) && getNaNPayload(x) == payload;
884 }
885
886 assert(feqrel(erfc(0.250L), erfc0_250 )>=real.mant_dig-1);
887 assert(feqrel(erfc(0.375L), erfc0_375 )>=real.mant_dig-0);
888 assert(feqrel(erfc(0.500L), erfc0_500 )>=real.mant_dig-2);
889 assert(feqrel(erfc(0.625L), erfc0_625 )>=real.mant_dig-1);
890 assert(feqrel(erfc(0.750L), erfc0_750 )>=real.mant_dig-1);
891 assert(feqrel(erfc(0.875L), erfc0_875 )>=real.mant_dig-4);
892 assert(feqrel(erfc(1.000L), erfc1_000 )>=real.mant_dig-2);
893 assert(feqrel(erfc(1.125L), erfc1_125 )>=real.mant_dig-2);
894 assert(feqrel(erf(0.875L), erf0_875 )>=real.mant_dig-1);
895 // The DMC implementation of erfc() fails this next test (just).
896 // Test point from Mathematica 11.0.
897 assert(feqrel(erfc(4.1L), 6.70002765408489837272673380763418472e-9L) >= real.mant_dig-5);
898
899 assert(isIdentical(erf(0.0),0.0));
900 assert(isIdentical(erf(-0.0),-0.0));
901 assert(erf(real.infinity) == 1.0);
902 assert(erf(-real.infinity) == -1.0);
903 assert(isNaNWithPayload(erf(NaN(0xDEF)), 0xDEF));
904 assert(isNaNWithPayload(erfc(NaN(0xDEF)), 0xDEF));
905 assert(isIdentical(erfc(real.infinity),0.0));
906 assert(erfc(-real.infinity) == 2.0);
907 assert(erfc(0) == 1.0);
908 }
909
910 /*
911 * Exponential of squared argument
912 *
913 * Computes y = exp(x*x) while suppressing error amplification
914 * that would ordinarily arise from the inexactness of the
915 * exponential argument x*x.
916 *
917 * If sign < 0, the result is inverted; i.e., y = exp(-x*x) .
918 *
919 * ACCURACY:
920 * Relative error:
921 * arithmetic domain # trials peak rms
922 * IEEE -106.566, 106.566 10^5 1.6e-19 4.4e-20
923 */
924
expx2(real x,int sign)925 real expx2(real x, int sign)
926 {
927 /*
928 Cephes Math Library Release 2.9: June, 2000
929 Copyright 2000 by Stephen L. Moshier
930 */
931 const real M = 32_768.0;
932 const real MINV = 3.0517578125e-5L;
933
934 x = abs(x);
935 if (sign < 0)
936 x = -x;
937
938 /* Represent x as an exact multiple of M plus a residual.
939 M is a power of 2 chosen so that exp(m * m) does not overflow
940 or underflow and so that |x - m| is small. */
941 real m = MINV * floor(M * x + 0.5L);
942 real f = x - m;
943
944 /* x^2 = m^2 + 2mf + f^2 */
945 real u = m * m;
946 real u1 = 2 * m * f + f * f;
947
948 if (sign < 0)
949 {
950 u = -u;
951 u1 = -u1;
952 }
953
954 if ((u+u1) > MAXLOG)
955 return real.infinity;
956
957 /* u is exact, u1 is small. */
958 return exp(u) * exp(u1);
959 }
960
961
962 /*
963 Computes the normal distribution function.
964
965 The normal (or Gaussian, or bell-shaped) distribution is
966 defined as:
967
968 normalDist(x) = 1/$(SQRT) π $(INTEGRAL -$(INFINITY), x) exp( - $(POWER t, 2)/2) dt
969 = 0.5 + 0.5 * erf(x/sqrt(2))
970 = 0.5 * erfc(- x/sqrt(2))
971
972 To maintain accuracy at high values of x, use
973 normalDistribution(x) = 1 - normalDistribution(-x).
974
975 Accuracy:
976 Within a few bits of machine resolution over the entire
977 range.
978
979 References:
980 $(LINK http://www.netlib.org/cephes/ldoubdoc.html),
981 G. Marsaglia, "Evaluating the Normal Distribution",
982 Journal of Statistical Software <b>11</b>, (July 2004).
983 */
normalDistributionImpl(real a)984 real normalDistributionImpl(real a)
985 {
986 real x = a * SQRT1_2;
987 real z = abs(x);
988
989 if ( z < 1.0 )
990 return 0.5L + 0.5L * erf(x);
991 else
992 {
993 real y = 0.5L * erfce(z);
994 /* Multiply by exp(-x^2 / 2) */
995 z = expx2(a, -1);
996 y = y * sqrt(z);
997 if ( x > 0.0L )
998 y = 1.0L - y;
999 return y;
1000 }
1001 }
1002
1003 @safe unittest
1004 {
1005 assert(fabs(normalDistributionImpl(1L) - (0.841344746068543))< 0.0000000000000005);
1006 assert(isIdentical(normalDistributionImpl(NaN(0x325)), NaN(0x325)));
1007 }
1008
1009 /*
1010 * Inverse of Normal distribution function
1011 *
1012 * Returns the argument, x, for which the area under the
1013 * Normal probability density function (integrated from
1014 * minus infinity to x) is equal to p.
1015 *
1016 * For small arguments 0 < p < exp(-2), the program computes
1017 * z = sqrt( -2 log(p) ); then the approximation is
1018 * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z) .
1019 * For larger arguments, x/sqrt(2 pi) = w + w^3 R(w^2)/S(w^2)) ,
1020 * where w = p - 0.5.
1021 */
1022 // TODO: isIEEEQuadruple (128 bit) real implementation; not available from CEPHES.
normalDistributionInvImpl(real p)1023 real normalDistributionInvImpl(real p)
1024 in {
1025 assert(p >= 0.0L && p <= 1.0L, "Domain error");
1026 }
1027 body
1028 {
1029 static immutable real[8] P0 =
1030 [ -0x1.758f4d969484bfdcp-7, 0x1.53cee17a59259dd2p-3,
1031 -0x1.ea01e4400a9427a2p-1, 0x1.61f7504a0105341ap+1, -0x1.09475a594d0399f6p+2,
1032 0x1.7c59e7a0df99e3e2p+1, -0x1.87a81da52edcdf14p-1, 0x1.1fb149fd3f83600cp-7
1033 ];
1034
1035 static immutable real[8] Q0 =
1036 [ -0x1.64b92ae791e64bb2p-7, 0x1.7585c7d597298286p-3,
1037 -0x1.40011be4f7591ce6p+0, 0x1.1fc067d8430a425ep+2, -0x1.21008ffb1e7ccdf2p+3,
1038 0x1.3d1581cf9bc12fccp+3, -0x1.53723a89fd8f083cp+2, 1.0
1039 ];
1040
1041 static immutable real[10] P1 =
1042 [ 0x1.20ceea49ea142f12p-13, 0x1.cbe8a7267aea80bp-7,
1043 0x1.79fea765aa787c48p-2, 0x1.d1f59faa1f4c4864p+1, 0x1.1c22e426a013bb96p+4,
1044 0x1.a8675a0c51ef3202p+5, 0x1.75782c4f83614164p+6, 0x1.7a2f3d90948f1666p+6,
1045 0x1.5cd116ee4c088c3ap+5, 0x1.1361e3eb6e3cc20ap+2
1046 ];
1047
1048 static immutable real[10] Q1 =
1049 [ 0x1.3a4ce1406cea98fap-13, 0x1.f45332623335cda2p-7,
1050 0x1.98f28bbd4b98db1p-2, 0x1.ec3b24f9c698091cp+1, 0x1.1cc56ecda7cf58e4p+4,
1051 0x1.92c6f7376bf8c058p+5, 0x1.4154c25aa47519b4p+6, 0x1.1b321d3b927849eap+6,
1052 0x1.403a5f5a4ce7b202p+4, 1.0
1053 ];
1054
1055 static immutable real[8] P2 =
1056 [ 0x1.8c124a850116a6d8p-21, 0x1.534abda3c2fb90bap-13,
1057 0x1.29a055ec93a4718cp-7, 0x1.6468e98aad6dd474p-3, 0x1.3dab2ef4c67a601cp+0,
1058 0x1.e1fb3a1e70c67464p+1, 0x1.b6cce8035ff57b02p+2, 0x1.9f4c9e749ff35f62p+1
1059 ];
1060
1061 static immutable real[8] Q2 =
1062 [ 0x1.af03f4fc0655e006p-21, 0x1.713192048d11fb2p-13,
1063 0x1.4357e5bbf5fef536p-7, 0x1.7fdac8749985d43cp-3, 0x1.4a080c813a2d8e84p+0,
1064 0x1.c3a4b423cdb41bdap+1, 0x1.8160694e24b5557ap+2, 1.0
1065 ];
1066
1067 static immutable real[8] P3 =
1068 [ -0x1.55da447ae3806168p-34, -0x1.145635641f8778a6p-24,
1069 -0x1.abf46d6b48040128p-17, -0x1.7da550945da790fcp-11, -0x1.aa0b2a31157775fap-8,
1070 0x1.b11d97522eed26bcp-3, 0x1.1106d22f9ae89238p+1, 0x1.029a358e1e630f64p+1
1071 ];
1072
1073 static immutable real[8] Q3 =
1074 [ -0x1.74022dd5523e6f84p-34, -0x1.2cb60d61e29ee836p-24,
1075 -0x1.d19e6ec03a85e556p-17, -0x1.9ea2a7b4422f6502p-11, -0x1.c54b1e852f107162p-8,
1076 0x1.e05268dd3c07989ep-3, 0x1.239c6aff14afbf82p+1, 1.0
1077 ];
1078
1079 if (p <= 0.0L || p >= 1.0L)
1080 {
1081 if (p == 0.0L)
1082 return -real.infinity;
1083 if ( p == 1.0L )
1084 return real.infinity;
1085 return real.nan; // domain error
1086 }
1087 int code = 1;
1088 real y = p;
1089 if ( y > (1.0L - EXP_2) )
1090 {
1091 y = 1.0L - y;
1092 code = 0;
1093 }
1094
1095 real x, z, y2, x0, x1;
1096
1097 if ( y > EXP_2 )
1098 {
1099 y = y - 0.5L;
1100 y2 = y * y;
1101 x = y + y * (y2 * rationalPoly( y2, P0, Q0));
1102 return x * SQRT2PI;
1103 }
1104
1105 x = sqrt( -2.0L * log(y) );
1106 x0 = x - log(x)/x;
1107 z = 1.0L/x;
1108 if ( x < 8.0L )
1109 {
1110 x1 = z * rationalPoly( z, P1, Q1);
1111 }
1112 else if ( x < 32.0L )
1113 {
1114 x1 = z * rationalPoly( z, P2, Q2);
1115 }
1116 else
1117 {
1118 x1 = z * rationalPoly( z, P3, Q3);
1119 }
1120 x = x0 - x1;
1121 if ( code != 0 )
1122 {
1123 x = -x;
1124 }
1125 return x;
1126 }
1127
1128
1129 @safe unittest
1130 {
1131 // TODO: Use verified test points.
1132 // The values below are from Excel 2003.
1133 assert(fabs(normalDistributionInvImpl(0.001) - (-3.09023230616779))< 0.00000000000005);
1134 assert(fabs(normalDistributionInvImpl(1e-50) - (-14.9333375347885))< 0.00000000000005);
1135 assert(feqrel(normalDistributionInvImpl(0.999), -normalDistributionInvImpl(0.001)) > real.mant_dig-6);
1136
1137 // Excel 2003 gets all the following values wrong!
1138 assert(normalDistributionInvImpl(0.0) == -real.infinity);
1139 assert(normalDistributionInvImpl(1.0) == real.infinity);
1140 assert(normalDistributionInvImpl(0.5) == 0);
1141 // (Excel 2003 returns norminv(p) = -30 for all p < 1e-200).
1142 // The value tested here is the one the function returned in Jan 2006.
1143 real unknown1 = normalDistributionInvImpl(1e-250L);
1144 assert( fabs(unknown1 -(-33.79958617269L) ) < 0.00000005);
1145 }
1146