@@ -1,4 +1,4 @@
1-/* statistics accelerator C extensor: _statistics module. */
1+/* statistics accelerator C extension: _statistics module. */
2233#include"Python.h"
44#include"structmember.h"
@@ -10,11 +10,13 @@ module _statistics
1010[clinic start generated code]*/
1111/*[clinic end generated code: output=da39a3ee5e6b4b0d input=864a6f59b76123b2]*/
121213-14-staticPyMethodDefspeedups_methods[] = {
15-_STATISTICS__NORMAL_DIST_INV_CDF_METHODDEF
16- {NULL, NULL, 0, NULL}
17-};
13+/*
14+ * There is no closed-form solution to the inverse CDF for the normal
15+ * distribution, so we use a rational approximation instead:
16+ * Wichura, M.J. (1988). "Algorithm AS241: The Percentage Points of the
17+ * Normal Distribution". Applied Statistics. Blackwell Publishing. 37
18+ * (3): 477–484. doi:10.2307/2347330. JSTOR 2347330.
19+ */
18201921/*[clinic input]
2022_statistics._normal_dist_inv_cdf -> double
@@ -34,7 +36,7 @@ _statistics__normal_dist_inv_cdf_impl(PyObject *module, double p, double mu,
3436// Algorithm AS 241: The Percentage Points of the Normal Distribution
3537if(fabs(q) <= 0.425) {
3638r=0.180625-q*q;
37-// Hash sum AB: 55.88319 28806 14901 4439
39+// Hash sum-55.8831928806149014439
3840num= (((((((2.5090809287301226727e+3*r+
39413.3430575583588128105e+4) *r+
40426.7265770927008700853e+4) *r+
@@ -54,11 +56,11 @@ _statistics__normal_dist_inv_cdf_impl(PyObject *module, double p, double mu,
5456x=num / den;
5557returnmu+ (x*sigma);
5658 }
57-r=q <= 0.0? p : 1.0-p;
59+r=(q <= 0.0) ? p : (1.0-p);
5860r=sqrt(-log(r));
5961if (r <= 5.0) {
6062r=r-1.6;
61-// Hash sum CD: 49.33206 50330 16102 89036
63+// Hash sum-49.33206503301610289036
6264num= (((((((7.74545014278341407640e-4*r+
63652.27238449892691845833e-2) *r+
64662.41780725177450611770e-1) *r+
@@ -77,7 +79,7 @@ _statistics__normal_dist_inv_cdf_impl(PyObject *module, double p, double mu,
77791.0);
7880 } else {
7981r-=5.0;
80-// Hash sum EF: 47.52583 31754 92896 71629
82+// Hash sum-47.52583317549289671629
8183num= (((((((2.01033439929228813265e-7*r+
82842.71155556874348757815e-5) *r+
83851.24266094738807843860e-3) *r+
@@ -96,23 +98,30 @@ _statistics__normal_dist_inv_cdf_impl(PyObject *module, double p, double mu,
96981.0);
9799 }
98100x=num / den;
99-if (q<0.0) x=-x;
101+if (q<0.0) {
102+x=-x;
103+ }
100104returnmu+ (x*sigma);
101105}
102106107+108+staticPyMethodDefstatistics_methods[] = {
109+_STATISTICS__NORMAL_DIST_INV_CDF_METHODDEF
110+ {NULL, NULL, 0, NULL}
111+};
112+103113staticstructPyModuleDefstatisticsmodule= {
104114PyModuleDef_HEAD_INIT,
105115"_statistics",
106116_statistics__normal_dist_inv_cdf__doc__,
107117-1,
108-speedups_methods,
118+statistics_methods,
109119NULL,
110120NULL,
111121NULL,
112122NULL
113123};
114124115-116125PyMODINIT_FUNC
117126PyInit__statistics(void)
118127{