diff --git a/scipy/_lib/boost_math b/scipy/_lib/boost_math
index 54e5acf7faea..260f7600e95d 160000
--- a/scipy/_lib/boost_math
+++ b/scipy/_lib/boost_math
@@ -1 +1 @@
-Subproject commit 54e5acf7faeab1400031d84bbfb2c46e81d82442
+Subproject commit 260f7600e95dcf87b228d46a4a37cf594df5437b
diff --git a/scipy/special/tests/test_basic.py b/scipy/special/tests/test_basic.py
index efb34b2c5238..42ee177fd630 100644
--- a/scipy/special/tests/test_basic.py
+++ b/scipy/special/tests/test_basic.py
@@ -1445,6 +1445,14 @@ def test_betainc_domain_errors(self, func, args):
             with pytest.raises(special.SpecialFunctionError, match='domain'):
                 special.betainc(*args)
 
+    @pytest.mark.parametrize('dtype', [np.float32, np.float64])
+    def test_gh21426(self, dtype):
+        # Test for gh-21426: betaincinv must not return NaN
+        a = np.array([5.], dtype=dtype)
+        x = np.array([0.5], dtype=dtype)
+        result = special.betaincinv(a, a, x)
+        assert_allclose(result, x, rtol=10 * np.finfo(dtype).eps)
+
 
 class TestCombinatorics:
     def test_comb(self):