/* statistics accelerator C extension: _statistics module. */

// Need limited C API version 3.13 for Py_mod_gil
#include "pyconfig.h"   // Py_GIL_DISABLED
#ifndef Py_GIL_DISABLED
#  define Py_LIMITED_API 0x030d0000
#endif

#include "Python.h"
#include "clinic/_statisticsmodule.c.h"

/*[clinic input]
module _statistics

[clinic start generated code]*/
/*[clinic end generated code: output=da39a3ee5e6b4b0d input=864a6f59b76123b2]*/

/*
 * There is no closed-form solution to the inverse CDF for the normal
 * distribution, so we use a rational approximation instead:
 * Wichura, M.J. (1988). "Algorithm AS241: The Percentage Points of the
 * Normal Distribution".  Applied Statistics. Blackwell Publishing. 37
 * (3): 477484. doi:10.2307/2347330. JSTOR 2347330.
 */

/*[clinic input]
_statistics._normal_dist_inv_cdf -> double
   p: double
   mu: double
   sigma: double
   /
[clinic start generated code]*/

static double
_statistics__normal_dist_inv_cdf_impl(PyObject *module, double p, double mu,
                                      double sigma)
/*[clinic end generated code: output=02fd19ddaab36602 input=24715a74be15296a]*/
{
    double q, num, den, r, x;
    if (p = 1.0) {
        goto error;
    }

    q = p - 0.5;
    if(fabs(q)