BUG: special: Fix unchecked memory allocations in specfun.h (scipy#22080)

* WIP: Fix unchecked memory allocation in aswfa()

Summary of the changes:

* Add another SF error code, "SF_ERROR_NOMEM", that is used by ufuncs
  that require internally allocated memory.  Just like the other possible
  errors, how this error is handled can be queried and set using `geterr()`
  and `seterr()`.  Unlike the other error conditions, the default behavior
  is to raise an exception if a memory allocation fails.

* Modify the function `aswfa()` in `specfun.h` to return a status code
  that indicates when a mempory allocation failed.

* Modify the functions `prolate_aswfa_nocv()`, `oblate_aswfa_nocv()`,
  `prolate_aswfa()` and `oblate_aswfa()`  in `sphd_wave.h` to check
  the return value of `aswfa()`, and set the error state and return values
  appropriately if `aswfa()` returns `NoMemory`.

* MAINT: Rename 'nomem' -> 'memory'.

* Fix all uses of unchecked allocations in specfun.h

* Fix aswfa() in specfunc.h

* Fix refguide check.

* Add comment about the return value when nonconvergence is detected in rmn2l.

* Change (c/m)alloc calls to use 'new' wrapped in a unique_ptr.

* Fix a few more names passed to set_error()

* Set outputs to nan in mtu12 when alloc fails
Also update docstring-like comments of some C++ functions.

Author: WarrenWeckesser

scipy/special/__init__.py

@@ -23,8 +23,9 @@
 
 Errors are handled by returning NaNs or other appropriate values.
 Some of the special function routines can emit warnings or raise
-exceptions when an error occurs. By default this is disabled; to
-query and control the current error handling state the following
+exceptions when an error occurs. By default this is disabled, except
+for memory allocation errors, which result in an exception being raised.
+To query and control the current error handling state the following
 functions are provided.
 
 .. autosummary::

scipy/special/_ellip_harm.pxd

@@ -98,7 +98,7 @@ cdef inline double* lame_coefficients(double h2, double k2, int n, int p,
                                + (sizeof(CBLAS_INT)*(2*size + liwork)))
     bufferp[0] = buffer
     if not buffer:
-        sf_error.error("ellip_harm", sf_error.NO_RESULT, "failed to allocate memory")
+        sf_error.error("ellip_harm", sf_error.MEMORY, "failed to allocate memory")
         return NULL
 
     cdef double *g = <double *>buffer
@@ -166,7 +166,7 @@ cdef inline double* lame_coefficients(double h2, double k2, int n, int p,
              &size, isuppz, work, &lwork, iwork, &liwork, &info)
 
     if info != 0:
-        sf_error.error("ellip_harm", sf_error.NO_RESULT, "failed to allocate memory")
+        sf_error.error("ellip_harm", sf_error.MEMORY, "failed to allocate memory")
         return NULL
 
     for i in range(0, size):

scipy/special/_specfun.pyx

@@ -16,22 +16,29 @@ cdef extern from "xsf/par_cyl.h" nogil:
     void specfun_pbdv 'xsf::detail::pbdv'(double x, double v, double *dv, double *dp, double *pdf, double *pdd)
     void specfun_pbvv 'xsf::detail::pbvv'(double x, double v, double *vv, double *vp, double *pvf, double *pvd)
 
+cdef extern from "xsf/specfun/specfun.h" namespace "xsf::specfun":
+
+    cpdef enum class Status:
+        OK = 0
+        NoMemory
+        Other
+
 cdef extern from "xsf/specfun/specfun.h" nogil:
     void specfun_bernob 'xsf::specfun::bernob'(int n, double *bn)
     void specfun_cerzo 'xsf::specfun::cerzo'(int nt, ccomplex[double] *zo)
     void specfun_cpbdn 'xsf::specfun::cpbdn'(int n, ccomplex[double] z, ccomplex[double] *cpb, ccomplex[double] *cpd)
     void specfun_cyzo 'xsf::specfun::cyzo'(int nt, int kf, int kc, ccomplex[double] *zo, ccomplex[double] *zv)
     void specfun_eulerb 'xsf::specfun::eulerb'(int n, double *en)
     void specfun_fcoef 'xsf::specfun::fcoef'(int kd, int m, double q, double a, double *fc)
-    void specfun_jdzo 'xsf::specfun::jdzo'(int nt, double *zo, int *n, int *m, int *p)
+    Status specfun_jdzo 'xsf::specfun::jdzo'(int nt, double *zo, int *n, int *m, int *p)
     void specfun_jyzo 'xsf::specfun::jyzo'(int n, int nt, double *rj0, double *rj1, double *ry0, double *ry1)
     void specfun_lamn 'xsf::specfun::lamn'(int n, double x, int *nm, double *bl, double *dl)
     void specfun_lamv 'xsf::specfun::lamv'(double v, double x, double *vm, double *vl, double *dl)
     void specfun_lqnb 'xsf::specfun::lqnb'(int n, double x, double* qn, double* qd)
     void specfun_pbdv 'xsf::specfun::pbdv'(double x, double v, double *dv, double *dp, double *pdf, double *pdd)
     void specfun_pbvv 'xsf::specfun::pbvv'(double x, double v, double *vv, double *vp, double *pvf, double *pvd)
-    void specfun_sdmn 'xsf::specfun::sdmn'(int m, int n, double c, double cv, double kd, double *df)
-    void specfun_segv 'xsf::specfun::segv'(int m, int n, double c, int kd, double *cv, double *eg)
+    Status specfun_sdmn 'xsf::specfun::sdmn'(int m, int n, double c, double cv, double kd, double *df)
+    Status specfun_segv 'xsf::specfun::segv'(int m, int n, double c, int kd, double *cv, double *eg)
 
 
 def airyzo(int nt, int kf):
@@ -224,7 +231,8 @@ def jdzo(int nt):
     mm = <int*>cnp.PyArray_DATA(m)
     pp = <int*>cnp.PyArray_DATA(p)
 
-    specfun_jdzo(nt, zzo, nn, mm, pp)
+    if (specfun_jdzo(nt, zzo, nn, mm, pp) == Status.NoMemory):
+        raise MemoryError('jnjnp_zeros: failed to allocate working memory in jdzo.')
 
     return n, m, p, zo
 
@@ -353,7 +361,8 @@ def sdmn(int m, int n, double c, double cv, int kd):
 
     df = cnp.PyArray_ZEROS(1, dims, cnp.NPY_FLOAT64, 0)
     ddf = <cnp.float64_t *>cnp.PyArray_DATA(df)
-    specfun_sdmn(m, n, c, cv, kd, ddf)
+    if specfun_sdmn(m, n, c, cv, kd, ddf) == Status.NoMemory:
+        raise MemoryError('sdmn: failed to allocate working memory.')
     return df
 
 
@@ -369,5 +378,11 @@ def segv(int m, int n, double c, int kd):
 
     eg = cnp.PyArray_ZEROS(1, dims, cnp.NPY_FLOAT64, 0)
     eeg = <cnp.float64_t *>cnp.PyArray_DATA(eg)
-    specfun_segv(m, n, c, kd, &cv, eeg)
+    if specfun_segv(m, n, c, kd, &cv, eeg) == Status.NoMemory:
+        # Note: segv is a private function that is called by either pro_cv_seq
+        # or obl_cv_seq.  We make the error message useful by including the
+        # approriate name in it.
+        caller = 'pro_cv_seq' if kd == 1 else 'obl_cv_seq'
+        msg = f'{caller}: failed to allocate working memory in segv.'
+        raise MemoryError(msg)
     return cv, eg

scipy/special/_special_ufuncs.cpp

@@ -802,7 +802,7 @@ PyMODINIT_FUNC PyInit__special_ufuncs() {
 
     PyObject *pro_cv = xsf::numpy::ufunc(
         {static_cast<xsf::numpy::fff_f>(xsf::prolate_segv), static_cast<xsf::numpy::ddd_d>(xsf::prolate_segv)},
-        "obl_cv", pro_cv_doc);
+        "pro_cv", pro_cv_doc);
     PyModule_AddObjectRef(_special_ufuncs, "pro_cv", pro_cv);
 
     PyObject *pro_rad1 = xsf::numpy::ufunc({static_cast<xsf::numpy::ffff_ff>(xsf::prolate_radial1_nocv),

scipy/special/_ufuncs_extra_code.pxi

@@ -12,7 +12,8 @@ _sf_error_code_map = {
     'no_result': 6,
     'domain': 7,
     'arg': 8,
-    'other': 9
+    'other': 9,
+    'memory': 10
 }
 
 _sf_error_action_map = {
@@ -54,11 +55,12 @@ def geterr():
 
     >>> import scipy.special as sc
     >>> for key, value in sorted(sc.geterr().items()):
-    ...     print("{}: {}".format(key, value))
+    ...     print(f'{key}: {value}')
     ...
     arg: ignore
     domain: ignore
     loss: ignore
+    memory: raise
     no_result: ignore
     other: ignore
     overflow: ignore

scipy/special/ellint_carlson_cpp_lite/ellint_common.hh

@@ -91,6 +91,7 @@ namespace ellint_carlson
 	bad_args,
 	bad_rerr,
 	other,
+	memory,
 	unused
     };
 

scipy/special/sf_error.cc

@@ -20,6 +20,7 @@ const char *sf_error_messages[] = {
     "domain error",
     "invalid input argument",
     "other error",
+    "memory allocation failed",
     NULL
 };
 
@@ -35,7 +36,7 @@ void sf_error_v(const char *func_name, sf_error_t code, const char *fmt, va_list
     static PyObject *py_SpecialFunctionWarning = NULL;
     sf_action_t action;
 
-    if ((int) code < 0 || (int) code >= 10) {
+    if ((int) code < 0 || (int) code >= SF_ERROR__LAST) {
         code = SF_ERROR_OTHER;
     }
     action = scipy_sf_error_get_action(code);

scipy/special/sf_error.pxd

@@ -12,6 +12,8 @@ cdef extern from "sf_error.h":
         DOMAIN "SF_ERROR_DOMAIN"
         ARG "SF_ERROR_ARG"
         OTHER "SF_ERROR_OTHER"
+        MEMORY "SF_ERROR_MEMORY"
+        LAST "SF_ERROR__LAST"
 
     ctypedef enum sf_action_t:
         IGNORE "SF_ERROR_IGNORE"
@@ -32,7 +34,7 @@ cdef inline int _sf_error_test_function(int code) noexcept nogil:
     """
     cdef sf_error_t sf_error_code
 
-    if code < 0 or code >= 10:
+    if code < 0 or code >= LAST:
         sf_error_code = OTHER
     else:
         sf_error_code = <sf_error_t>code

scipy/special/sf_error_state.c

@@ -15,6 +15,7 @@ static volatile sf_action_t sf_error_actions[] = {
     SF_ERROR_IGNORE, /* SF_ERROR_DOMAIN */
     SF_ERROR_IGNORE, /* SF_ERROR_ARG */
     SF_ERROR_IGNORE, /* SF_ERROR_OTHER */
+    SF_ERROR_RAISE,  /* SF_ERROR_MEMORY */
     SF_ERROR_IGNORE  /* SF_ERROR__LAST */
 };
 

scipy/special/tests/test_sf_error.py

@@ -19,7 +19,8 @@
     'no_result': 6,
     'domain': 7,
     'arg': 8,
-    'other': 9
+    'other': 9,
+    'memory': 10,
 }
 
 _sf_error_actions = [