Repository: scipy/scipy
Pull request: BUG: special: Use Boost for ncfdtr, fixing accuracy issues
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
https://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message.
However, please only include an issue number in the description, not the title,
and please ensure that any code names containing underscores are enclosed in backticks.

Depending on your changes, you can skip CI operations and save time and energy: 
https://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Closes #15582, ticks a box for #20223 

#### What does this implement/fix?
<!--Please explain your changes.-->
This PR changes the scalar kernel for `ncfdtr` to use a wrapper of the more accurate Boost implementation which was already being used in stats for `ncf.cdf`. Previously there were two ufuncs a public `ncfdtr` using cdflib, and a private `_ncf_cdf` using Boost. I've changed it so there is no only one ufunc, `ncfdtr` which uses Boost. I had to change the argument order for `_ncf_cdf`, which differed from that of `ncfdtr`, but this was not difficult because `_ncf_cdf` is only used in one place in stats. I also had to add `NAN` handling in the Boost wrapper.

I had to `xfail` round trip tests for `ncfdtri`, `ncfdtridfn`, `ncfdtridfn`, and `ncfdtrinc`, because these functions invert the less accurate `cdflib` implementation of `ncfdtr`, and therefore now fail. We can fix this after https://github.com/scipy/scipy/pull/21454 lands and we are able to have these functions invert the Boost implementation now being used.

I added some tests comparing `ncfdtr` against an `mpmath` reference. 11 of the 12 test cases I added fail on main. Locally I tested a larger range consisting of 625 cases. On `main` the deciles for relative error over these cases are

```
array([0.00000000e+00, 1.15300845e-14, 4.37093510e-13, 5.72470839e-10,
       1.48813582e-08, 7.50889455e-08, 2.50213950e-07, 1.19707132e-06,
       1.48768369e-05, 3.80639979e-05, 2.12464729e-04])
```

using Boost, the deciles for relative error are

```
array([0.00000000e+00, 0.00000000e+00, 1.40500993e-16, 2.52732515e-16,
       4.34020841e-16, 6.33097261e-16, 1.09243091e-15, 1.53843495e-15,
       2.34037246e-15, 6.13809455e-15, 8.41339873e-14])
```

The max relative error for Boost is better than the 20th percentile for cdflib.

The script generating the test cases is in the comments for the test I added to `special/tests/test_cdflib.py`.
#### Additional information
As suggested by @dschmitz89 here, https://github.com/scipy/scipy/issues/15582#issuecomment-2177581310, we should follow up by using the Boost implementations instead of cdflib when stats is already using Boost, and investigate Boost or other alternatives to cdflib for other functions.)

<!--Any additional information you think is important.-->


Problem statement:
BUG: `special.ncfdtr(dfn, dfd, nc, f)`  and `stats.ncf.cdf(x, dfn, dfd, nc)` results are different from R and Matlab results.
### Describe your issue.


Calculated results are not equal to that of R and Matlab

`stats.ncf.cdf(20, 6,33,30.4) = '
0.998868`

but in R
`pf(20, df1 = 6, df2 = 33, ncp = 30.4) = 
0.998921`

In Matlab
`ncfcdf(20,  6, 33, 30.4) =
0.998921`


### Reproducing Code Example

```python
No
```


### Error message

```shell
No
```


### SciPy/NumPy/Python version information

scipy.__version__ '1.6.2'

Hints from issue discussion:
SciPy 1.8.0 outputs same value.
```
In [7]: import sys, scipy, numpy;print(scipy.__version__, numpy.__version__, sys.version_info)
1.8.0 1.22.2 sys.version_info(major=3, minor=9, micro=2, releaselevel='final', serial=0)

In [8]: scipy.special.ncfdtr(6, 33, 30.4, 20)
Out[8]: 0.9988687760741152

In [9]: scipy.stats.ncf.cdf(20, 6, 33, 30.4)
Out[9]: 0.9988687760741152
```
Another one to consider for Boost, @mckib2?
Actually we only fixed this for `stats`, though. It should probably be fixed for special.
Good idea. Boost is very good for non-central distributions (for example R gives wrong values for the non-central Student when nc > 37.62).
To resolve this and reduce the number of duplications across scipy, could we just let `ncfdtr` (and friends such as `nctdtr`, `nbdtr` ...) dispatch to the boost functions? I have a strong opinion that we should remove CDFlib functions when better alternatives exist.
Is this still an issue ? I get the following on main

```python
In [2]: stats.ncf.cdf(20, 6, 33, 30.4)
Out[2]: np.float64(0.9989211030878699)
```
Ah OK special is still broken nevermind the noise

Commit message: Fix references to stats dist methods in docstrings [docs only]

Patch to review:
diff --git a/scipy/special/_add_newdocs.py b/scipy/special/_add_newdocs.py
index 409eb64b4885..2dc3d9d63cba 100644
--- a/scipy/special/_add_newdocs.py
+++ b/scipy/special/_add_newdocs.py
@@ -6821,7 +6821,7 @@ def add_newdoc(name, doc):
     dfd : array_like
         Degrees of freedom of the denominator sum of squares.  Range (0, inf).
     nc : array_like
-        Noncentrality parameter.  Should be in range (0, 1e4).
+        Noncentrality parameter.  Range [0, inf).
     f : array_like
         Quantiles, i.e. the upper limit of integration.
     out : ndarray, optional
@@ -6839,10 +6839,12 @@ def add_newdoc(name, doc):
     ncfdtridfd : Inverse of `ncfdtr` with respect to `dfd`.
     ncfdtridfn : Inverse of `ncfdtr` with respect to `dfn`.
     ncfdtrinc : Inverse of `ncfdtr` with respect to `nc`.
+    scipy.stats.ncf : Non-central F distribution.
 
     Notes
     -----
-    Wrapper for the CDFLIB [1]_ Fortran routine `cdffnc`.
+    This function calculates the CDF of the non-central f distribution using
+    the Boost Math C++ library [1]_.
 
     The cumulative distribution function is computed using Formula 26.6.20 of
     [2]_:
@@ -6854,16 +6856,13 @@ def add_newdoc(name, doc):
     where :math:`I` is the regularized incomplete beta function, and
     :math:`x = f d_n/(f d_n + d_d)`.
 
-    The computation time required for this routine is proportional to the
-    noncentrality parameter `nc`.  Very large values of this parameter can
-    consume immense computer resources.  This is why the search range is
-    bounded by 10,000.
+    Note that argument order of `ncfdtr` is different from that of the
+    otherwise equivalent ``cdf`` method of `scipy.stats.ncf`. `f` is the last
+    parameter of `ncfdtr` but the first parameter of ``scipy.stats.ncf.cdf``.
 
     References
     ----------
-    .. [1] Barry Brown, James Lovato, and Kathy Russell,
-           CDFLIB: Library of Fortran Routines for Cumulative Distribution
-           Functions, Inverses, and Other Parameters.
+    .. [1] The Boost Developers. "Boost C++ Libraries". https://www.boost.org/.
     .. [2] Milton Abramowitz and Irene A. Stegun, eds.
            Handbook of Mathematical Functions with Formulas,
            Graphs, and Mathematical Tables. New York: Dover, 1972.
@@ -6907,7 +6906,7 @@ def add_newdoc(name, doc):
     dfd : array_like
         Degrees of freedom of the denominator sum of squares.  Range (0, inf).
     nc : array_like
-        Noncentrality parameter.  Should be in range (0, 1e4).
+        Noncentrality parameter.  Range [0, inf).
     p : array_like
         Value of the cumulative distribution function.  Must be in the
         range [0, 1].
@@ -6925,6 +6924,20 @@ def add_newdoc(name, doc):
     ncfdtridfd : Inverse of `ncfdtr` with respect to `dfd`.
     ncfdtridfn : Inverse of `ncfdtr` with respect to `dfn`.
     ncfdtrinc : Inverse of `ncfdtr` with respect to `nc`.
+    scipy.stats.ncf : Non-central F distribution.
+
+    Notes
+    -----
+    This function calculates the Quantile of the non-central f distribution
+    using the Boost Math C++ library [1]_.
+
+    Note that argument order of `ncfdtri` is different from that of the
+    otherwise equivalent ``ppf`` method of `scipy.stats.ncf`. `p` is the last
+    parameter of `ncfdtri` but the first parameter of ``scipy.stats.ncf.ppf``.
+
+    References
+    ----------
+    .. [1] The Boost Developers. "Boost C++ Libraries". https://www.boost.org/.
 
     Examples
     --------
diff --git a/scipy/special/boost_special_functions.h b/scipy/special/boost_special_functions.h
index 65f6d12b80a5..f58f0fbddcd0 100644
--- a/scipy/special/boost_special_functions.h
+++ b/scipy/special/boost_special_functions.h
@@ -760,46 +760,86 @@ ncf_pdf_double(double x, double v1, double v2, double l)
 
 template<typename Real>
 Real
-ncf_cdf_wrap(const Real x, const Real v1, const Real v2, const Real l)
+ncf_cdf_wrap(const Real v1, const Real v2, const Real l, const Real x)
 {
-    if (std::isfinite(x)) {
-        return boost::math::cdf(
-            boost::math::non_central_f_distribution<Real, StatsPolicy>(v1, v2, l), x);
+    if (std::isnan(x) || std::isnan(v1) || std::isnan(v2) || std::isnan(l)) {
+	return NAN;
     }
-    // -inf => 0, inf => 1
-    return 1.0 - std::signbit(x);
+    if ((v1 <= 0) || (v2 <= 0) || (l < 0)) {
+	sf_error("ncfdtr", SF_ERROR_DOMAIN, NULL);
+	return NAN;
+    }
+    if (std::isinf(x)) {
+	// -inf => 0, inf => 1
+	return 1.0 - std::signbit(x);
+    }
+    Real y;
+    try {
+	y = boost::math::cdf(
+                boost::math::non_central_f_distribution<Real, SpecialPolicy>(v1, v2, l), x);
+    } catch (...) {
+	/* Boost was unable to produce a result. Can happen when one or both
+	 * of v1 and v2 is very small and x is very large. e.g.
+	 * ncdtr(1e-100, 3, 1.5, 1e100). */
+        sf_error("ncfdtr", SF_ERROR_NO_RESULT, NULL);
+        y = NAN;
+    }
+    return y;
 }
 
 float
-ncf_cdf_float(float x, float v1, float v2, float l)
+ncf_cdf_float(float v1, float v2, float l, float x)
 {
-    return ncf_cdf_wrap(x, v1, v2, l);
+    return ncf_cdf_wrap(v1, v2, l, x);
 }
 
 double
-ncf_cdf_double(double x, double v1, double v2, double l)
+ncf_cdf_double(double v1, double v2, double l, double x)
 {
-    return ncf_cdf_wrap(x, v1, v2, l);
+    return ncf_cdf_wrap(v1, v2, l, x);
 }
 
 template<typename Real>
 Real
-ncf_ppf_wrap(const Real x, const Real v1, const Real v2, const Real l)
+ncf_ppf_wrap(const Real v1, const Real v2, const Real l, const Real x)
 {
-    return boost::math::quantile<Real, StatsPolicy>(
-        boost::math::non_central_f_distribution<Real, StatsPolicy>(v1, v2, l), x);
+    if (std::isnan(x) || std::isnan(v1) || std::isnan(v2) || std::isnan(l)) {
+	return NAN;
+    }
+    if ((v1 <= 0) || (v2 <= 0) || (l < 0) || (x < 0) || (x > 1)) {
+	sf_error("ncfdtr", SF_ERROR_DOMAIN, NULL);
+	return NAN;
+    }
+    Real y;
+    try {
+	y = boost::math::quantile<Real, SpecialPolicy>(
+                boost::math::non_central_f_distribution<Real, SpecialPolicy>(v1, v2, l), x);
+    } catch (const std::domain_error& e) {
+        sf_error("ncfdtri", SF_ERROR_DOMAIN, NULL);
+        y = NAN;
+    } catch (const std::overflow_error& e) {
+        sf_error("ncfdtri", SF_ERROR_OVERFLOW, NULL);
+        y = INFINITY;
+    } catch (const std::underflow_error& e) {
+        sf_error("ncfdtri", SF_ERROR_UNDERFLOW, NULL);
+        y = 0;
+    } catch (...) {
+        sf_error("ncfdtri", SF_ERROR_NO_RESULT, NULL);
+        y = NAN;
+    }
+    return y;
 }
 
 float
-ncf_ppf_float(float x, float v1, float v2, float l)
+ncf_ppf_float(float v1, float v2, float l, float x)
 {
-    return ncf_ppf_wrap(x, v1, v2, l);
+    return ncf_ppf_wrap(v1, v2, l, x);
 }
 
 double
-ncf_ppf_double(double x, double v1, double v2, double l)
+ncf_ppf_double(double v1, double v2, double l, double x)
 {
-    return ncf_ppf_wrap(x, v1, v2, l);
+    return ncf_ppf_wrap(v1, v2, l, x);
 }
 
 template<typename Real>
diff --git a/scipy/special/cython_special.pxd b/scipy/special/cython_special.pxd
index a472206d75fd..add0329ac05b 100644
--- a/scipy/special/cython_special.pxd
+++ b/scipy/special/cython_special.pxd
@@ -190,8 +190,8 @@ cpdef double nbdtrc(dlp_number_t x0, dlp_number_t x1, double x2) noexcept nogil
 cpdef double nbdtri(dlp_number_t x0, dlp_number_t x1, double x2) noexcept nogil
 cpdef double nbdtrik(double x0, double x1, double x2) noexcept nogil
 cpdef double nbdtrin(double x0, double x1, double x2) noexcept nogil
-cpdef double ncfdtr(double x0, double x1, double x2, double x3) noexcept nogil
-cpdef double ncfdtri(double x0, double x1, double x2, double x3) noexcept nogil
+cpdef df_number_t ncfdtr(df_number_t x0, df_number_t x1, df_number_t x2, df_number_t x3) noexcept nogil
+cpdef df_number_t ncfdtri(df_number_t x0, df_number_t x1, df_number_t x2, df_number_t x3) noexcept nogil
 cpdef double ncfdtridfd(double x0, double x1, double x2, double x3) noexcept nogil
 cpdef double ncfdtridfn(double x0, double x1, double x2, double x3) noexcept nogil
 cpdef double ncfdtrinc(double x0, double x1, double x2, double x3) noexcept nogil
@@ -258,4 +258,4 @@ cpdef Dd_number_t yve(double x0, Dd_number_t x1) noexcept nogil
 cpdef double zetac(double x0) noexcept nogil
 cpdef double wright_bessel(double x0, double x1, double x2) noexcept nogil
 cpdef double log_wright_bessel(double x0, double x1, double x2) noexcept nogil
-cpdef double ndtri_exp(double x0) noexcept nogil
\ No newline at end of file
+cpdef double ndtri_exp(double x0) noexcept nogil
diff --git a/scipy/special/cython_special.pyx b/scipy/special/cython_special.pyx
index 28076bc08c95..6681f26884f7 100644
--- a/scipy/special/cython_special.pyx
+++ b/scipy/special/cython_special.pyx
@@ -1710,14 +1710,6 @@ from ._cdflib_wrappers cimport nbdtrin as _func_nbdtrin
 ctypedef double _proto_nbdtrin_t(double, double, double) noexcept nogil
 cdef _proto_nbdtrin_t *_proto_nbdtrin_t_var = &_func_nbdtrin
 
-from ._cdflib_wrappers cimport ncfdtr as _func_ncfdtr
-ctypedef double _proto_ncfdtr_t(double, double, double, double) noexcept nogil
-cdef _proto_ncfdtr_t *_proto_ncfdtr_t_var = &_func_ncfdtr
-
-from ._cdflib_wrappers cimport ncfdtri as _func_ncfdtri
-ctypedef double _proto_ncfdtri_t(double, double, double, double) noexcept nogil
-cdef _proto_ncfdtri_t *_proto_ncfdtri_t_var = &_func_ncfdtri
-
 from ._cdflib_wrappers cimport ncfdtridfd as _func_ncfdtridfd
 ctypedef double _proto_ncfdtridfd_t(double, double, double, double) noexcept nogil
 cdef _proto_ncfdtridfd_t *_proto_ncfdtridfd_t_var = &_func_ncfdtridfd
@@ -3175,13 +3167,29 @@ cpdef double nbdtrin(double x0, double x1, double x2) noexcept nogil:
     """See the documentation for scipy.special.nbdtrin"""
     return _func_nbdtrin(x0, x1, x2)
 
-cpdef double ncfdtr(double x0, double x1, double x2, double x3) noexcept nogil:
+cpdef df_number_t ncfdtr(df_number_t x0, df_number_t x1, df_number_t x2, df_number_t x3) noexcept nogil:
     """See the documentation for scipy.special.ncfdtr"""
-    return _func_ncfdtr(x0, x1, x2, x3)
+    if df_number_t is float:
+        return (<float(*)(float, float, float, float) noexcept nogil>scipy.special._ufuncs_cxx._export_ncf_cdf_float)(x0, x1, x2, x3)
+    elif df_number_t is double:
+        return (<double(*)(double, double, double, double) noexcept nogil>scipy.special._ufuncs_cxx._export_ncf_cdf_double)(x0, x1, x2, x3)
+    else:
+        if df_number_t is double:
+            return NAN
+        else:
+            return NAN
 
-cpdef double ncfdtri(double x0, double x1, double x2, double x3) noexcept nogil:
+cpdef df_number_t ncfdtri(df_number_t x0, df_number_t x1, df_number_t x2, df_number_t x3) noexcept nogil:
     """See the documentation for scipy.special.ncfdtri"""
-    return _func_ncfdtri(x0, x1, x2, x3)
+    if df_number_t is float:
+        return (<float(*)(float, float, float, float) noexcept nogil>scipy.special._ufuncs_cxx._export_ncf_ppf_float)(x0, x1, x2, x3)
+    elif df_number_t is double:
+        return (<double(*)(double, double, double, double) noexcept nogil>scipy.special._ufuncs_cxx._export_ncf_ppf_double)(x0, x1, x2, x3)
+    else:
+        if df_number_t is double:
+            return NAN
+        else:
+            return NAN
 
 cpdef double ncfdtridfd(double x0, double x1, double x2, double x3) noexcept nogil:
     """See the documentation for scipy.special.ncfdtridfd"""
diff --git a/scipy/special/functions.json b/scipy/special/functions.json
index 0b75ceb002f5..a5c1f1272f80 100644
--- a/scipy/special/functions.json
+++ b/scipy/special/functions.json
@@ -598,13 +598,15 @@
         }
     },
     "ncfdtr": {
-        "_cdflib_wrappers.pxd": {
-            "ncfdtr": "dddd->d"
+        "boost_special_functions.h++": {
+            "ncf_cdf_float": "ffff->f",
+            "ncf_cdf_double": "dddd->d"
         }
     },
     "ncfdtri": {
-        "_cdflib_wrappers.pxd": {
-            "ncfdtri": "dddd->d"
+        "boost_special_functions.h++": {
+            "ncf_ppf_float": "ffff->f",
+            "ncf_ppf_double": "dddd->d"
         }
     },
     "ncfdtridfd": {
@@ -910,18 +912,6 @@
             "ncf_pdf_double": "dddd->d"
         }
     },
-    "_ncf_cdf": {
-        "boost_special_functions.h++": {
-            "ncf_cdf_float": "ffff->f",
-            "ncf_cdf_double": "dddd->d"
-        }
-    },
-    "_ncf_ppf": {
-        "boost_special_functions.h++": {
-            "ncf_ppf_float": "ffff->f",
-            "ncf_ppf_double": "dddd->d"
-        }
-    },
     "_ncf_sf": {
         "boost_special_functions.h++": {
             "ncf_sf_float": "ffff->f",
diff --git a/scipy/special/tests/test_basic.py b/scipy/special/tests/test_basic.py
index cb148b09d3a6..efb34b2c5238 100644
--- a/scipy/special/tests/test_basic.py
+++ b/scipy/special/tests/test_basic.py
@@ -722,16 +722,34 @@ def test_ncfdtri(self):
         p = cephes.ncfdtr(2, 3, 1.5, f)
         assert_allclose(cephes.ncfdtri(2, 3, 1.5, p), f)
 
+    @pytest.mark.xfail(
+        reason=(
+            "ncfdtr uses a Boost math implementation but ncfdtridfd"
+            "inverts the less accurate cdflib implementation of ncfdtr."
+        )
+    )
     def test_ncfdtridfd(self):
         dfd = [1, 2, 3]
         p = cephes.ncfdtr(2, dfd, 0.25, 15)
         assert_allclose(cephes.ncfdtridfd(2, p, 0.25, 15), dfd)
 
+    @pytest.mark.xfail(
+        reason=(
+            "ncfdtr uses a Boost math implementation but ncfdtridfn"
+            "inverts the less accurate cdflib implementation of ncfdtr."
+        )
+    )
     def test_ncfdtridfn(self):
         dfn = [0.1, 1, 2, 3, 1e4]
         p = cephes.ncfdtr(dfn, 2, 0.25, 15)
         assert_allclose(cephes.ncfdtridfn(p, 2, 0.25, 15), dfn, rtol=1e-5)
 
+    @pytest.mark.xfail(
+        reason=(
+            "ncfdtr uses a Boost math implementation but ncfdtrinc"
+            "inverts the less accurate cdflib implementation of ncfdtr."
+        )
+    )
     def test_ncfdtrinc(self):
         nc = [0.5, 1.5, 2.0]
         p = cephes.ncfdtr(2, 3, nc, 15)
diff --git a/scipy/special/tests/test_cdflib.py b/scipy/special/tests/test_cdflib.py
index ca3e82299824..08e34f158de4 100644
--- a/scipy/special/tests/test_cdflib.py
+++ b/scipy/special/tests/test_cdflib.py
@@ -3,7 +3,6 @@
 
 The following functions still need tests:
 
-- ncfdtr
 - ncfdtri
 - ncfdtridfn
 - ncfdtridfd
@@ -525,3 +524,64 @@ def test_bdtrik_nbdtrik_inf():
         [np.nan, -np.inf, -10.0, -1.0, 0.0, .00001, .5, 1.0, np.inf])
     assert np.all(np.isnan(sp.bdtrik(y, np.inf, p)))
     assert np.all(np.isnan(sp.nbdtrik(y, np.inf, p)))
+
+
+@pytest.mark.parametrize(
+    "dfn,dfd,nc,f,expected",
+    [[100.0, 0.1, 0.1, 100.0, 0.29787396410092676],
+     [100.0, 100.0, 0.01, 0.1, 4.4344737598690424e-26],
+     [100.0, 0.01, 0.1, 0.01, 0.002848616633080384],
+     [10.0, 0.01, 1.0, 0.1, 0.012339557729057956],
+     [100.0, 100.0, 0.01, 0.01, 1.8926477420964936e-72],
+     [1.0, 100.0, 100.0, 0.1, 1.7925940526821304e-22],
+     [1.0, 0.01, 100.0, 10.0, 0.012334711965024968],
+     [1.0, 0.01, 10.0, 0.01, 0.00021944525290299],
+     [10.0, 1.0, 0.1, 100.0, 0.9219345555070705],
+     [0.1, 0.1, 1.0, 1.0, 0.3136335813423239],
+     [100.0, 100.0, 0.1, 10.0, 1.0],
+     [1.0, 0.1, 100.0, 10.0, 0.02926064279680897]]
+)
+def test_ncfdtr(dfn, dfd, nc, f, expected):
+    # Reference values computed with mpmath with the following script
+    #
+    # import numpy as np
+    #
+    # from mpmath import mp
+    # from scipy.special import ncfdtr
+    #
+    # mp.dps = 100
+    #
+    # def mp_ncfdtr(dfn, dfd, nc, f):
+    #     # Uses formula 26.2.20 from Abramowitz and Stegun.
+    #     dfn, dfd, nc, f = map(mp.mpf, (dfn, dfd, nc, f))
+    #     def term(j):
+    #         result = mp.exp(-nc/2)*(nc/2)**j / mp.factorial(j)
+    #         result *= mp.betainc(
+    #             dfn/2 + j, dfd/2, 0, f*dfn/(f*dfn + dfd), regularized=True
+    #         )
+    #         return result
+    #     result = mp.nsum(term, [0, mp.inf])
+    #     return float(result)
+    #
+    # dfn = np.logspace(-2, 2, 5)
+    # dfd = np.logspace(-2, 2, 5)
+    # nc = np.logspace(-2, 2, 5)
+    # f = np.logspace(-2, 2, 5)
+    #
+    # dfn, dfd, nc, f = np.meshgrid(dfn, dfd, nc, f)
+    # dfn, dfd, nc, f = map(np.ravel, (dfn, dfd, nc, f))
+    #
+    # cases = []
+    # re = []
+    # for x0, x1, x2, x3 in zip(*(dfn, dfd, nc, f)):
+    #     observed = ncfdtr(x0, x1, x2, x3)
+    #     expected = mp_ncfdtr(x0, x1, x2, x3)
+    #     cases.append((x0, x1, x2, x3, expected))
+    #     re.append((abs(expected - observed)/abs(expected)))
+    #
+    # assert np.max(re) < 1e-13
+    #
+    # rng = np.random.default_rng(1234)
+    # sample_idx = rng.choice(len(re), replace=False, size=12)
+    # cases = np.array(cases)[sample_idx].tolist()
+    assert_allclose(sp.ncfdtr(dfn, dfd, nc, f), expected, rtol=1e-13, atol=0)
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index e12fb23df27f..52ea41a3bd94 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -7718,11 +7718,11 @@ def _pdf(self, x, dfn, dfd, nc):
         return scu._ncf_pdf(x, dfn, dfd, nc)
 
     def _cdf(self, x, dfn, dfd, nc):
-        return scu._ncf_cdf(x, dfn, dfd, nc)
+        return sc.ncfdtr(dfn, dfd, nc, x)
 
     def _ppf(self, q, dfn, dfd, nc):
         with np.errstate(over='ignore'):  # see gh-17432
-            return scu._ncf_ppf(q, dfn, dfd, nc)
+            return sc.ncfdtri(dfn, dfd, nc, q)
 
     def _sf(self, x, dfn, dfd, nc):
         return scu._ncf_sf(x, dfn, dfd, nc)
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index b9207cfc6fcd..287d6225e9e3 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -9157,16 +9157,13 @@ def test_ncf_cdf_spotcheck():
     assert_allclose(check_val, np.round(scipy_val, decimals=6))
 
 
-@pytest.mark.skipif(sys.maxsize <= 2**32,
-                    reason="On some 32-bit the warning is not raised")
 def test_ncf_ppf_issue_17026():
     # Regression test for gh-17026
     x = np.linspace(0, 1, 600)
     x[0] = 1e-16
     par = (0.1, 2, 5, 0, 1)
-    with pytest.warns(RuntimeWarning):
-        q = stats.ncf.ppf(x, *par)
-        q0 = [stats.ncf.ppf(xi, *par) for xi in x]
+    q = stats.ncf.ppf(x, *par)
+    q0 = [stats.ncf.ppf(xi, *par) for xi in x]
     assert_allclose(q, q0)
 
 


Code review comments:
diff_hunk: @@ -760,46 +760,86 @@ ncf_pdf_double(double x, double v1, double v2, double l)
 
 template<typename Real>
 Real
-ncf_cdf_wrap(const Real x, const Real v1, const Real v2, const Real l)
+ncf_cdf_wrap(const Real v1, const Real v2, const Real l, const Real x)
 {
-    if (std::isfinite(x)) {
-        return boost::math::cdf(
-            boost::math::non_central_f_distribution<Real, StatsPolicy>(v1, v2, l), x);
+    if (std::isnan(x) || std::isnan(v1) || std::isnan(v2) || std::isnan(l)) {
+	return NAN;
     }
-    // -inf => 0, inf => 1
-    return 1.0 - std::signbit(x);
+    if ((v1 <= 0) || (v2 <= 0) || (l < 0)) {
+	sf_error("ncfdtr", SF_ERROR_DOMAIN, NULL);
+	return NAN;
+    }
+    if (std::isinf(x)) {
+	// -inf => 0, inf => 1
path: scipy/special/boost_special_functions.h
text: @user1:
To be consistent with previous behavior, this should return NaN for all `x < 0`, right?

That said, I was surprised when I realized that is how the special functions that are supposed to be distribution functions behave. But I think consistency is important, so I'd rather not change *only* the `ncf` functions 

@author:
Ah, good point. That trips me up too. I'll change it. I think I understand the rationale at least. For example, if we defined `betainc` to be the CDF of the Beta distribution, then it's natural to say `betainc(2, 2, -0.2)` is zero, since if $X \sim B(2, 2)$, then $P(X < -0.2) = 0$. Treating it just as a special function though, the natural way to extend things is by analytic continuation, and if we haven't implemented the analytic continuation (or if we're looking at the real valued version and the analytic continuation on the real axis but outside of the usual domain for the distribution function returns complex values) then it makes sense to return `NAN`.

@user1:
> For example, if we defined betainc to be the CDF of the Beta distribution

I think that's my point, though: `betainc` *isn't*/wasn't defined as the CDF of the beta distribution; `btdtr` is / was. So I would extend `betainc`'s domain (or returned NaN), whereas I would think `btdtr` would only be for reals and would return 0 or 1 outside the distribution's domain.

@author:
Yeah, I agree it's weird that we have different versions of the some functions, one a CDF, the other not, but the CDF doesn't behave as one would expect. I would imagine it's just an historical thing. It seems like the consistent principle being followed may be "whether something says it's a CDF or not, we don't treat it as a CDF, and the only valid way to extend a domain is by analytic continuation"?

Would you support maintaining CDF and non-CDF versions of the same functions with the only difference being how they treat arguments outside of the usual domain? It seems like a reasonable thing to do.

@user1:
> It seems like the consistent principle being followed may be "whether something says it's a CDF or not, we don't treat it as a CDF, and the only valid way to extend a domain is by analytic continuation"?

I think that describes what is happening. And with that in mind, yes, I imagine the reason for the apparent duplication is just historical.

> Would you support maintaining CDF and non-CDF versions of the same functions with the only difference being how they treat arguments outside of the usual domain?

Personally, I wouldn't want to work toward it because I think the random variable infrastructure can serve the need for functions to be treated as CDFs. On one hand, it would be great if the infrastructure could just call the special function and it would behave as a CDF, etc. There is lower overhead if this behavior happens at the compiled level. But if there is a need for this sort of behavior to be added by the distribution infrastructure (e.g. user-defined distributions) it might as well be applied consistently to all distributions, and then there is not much need to duplicate it in `special`. This also helps keep the usage distinction between `special` and `stats` obvious.

That said, I'm not sure what is best when there is already duplication or when there is a CDF for which analytic continuation is not useful.