fix: prevent array out-of-bounds write in segv function.Add boundary …

[Inuyasha-fan]

For details, please see: scipy/scipy#23973 .

segv Memory Corruption Analysis

1. Root Cause of Memory Corruption

There is an array out-of-bounds write in the segv function in specfun.h source code at lines 6058-6061.

// specfun.h segv
template <typename T> Status segv(int m, int n, T c, int kd, T *cv, T *eg) {
    // ...
    // 6058 - 6061
    if (l == 0)
        eg[2 * k - 2] = cv0[k - 1];
    if (l == 1)
        eg[2 * k - 1] = cv0[k - 1];
    // ...
}

The length of the eg array is n - m + 1, which comes from the segv function call in _specfun.pyx to the C++ function. In the segv function in specfun.h, the index of eg can exactly reach n - m + 1, causing an out-of-bounds write.

def segv(int m, int n, double c, int kd):
    """
    Compute the characteristic values of spheroidal wave functions.
    This is a wrapper for the function 'specfun_segv'.
    """
    cdef double cv
    cdef double *eeg
    cdef cnp.npy_intp dims[1]
    dims[0] = n - m + 1

    eg = cnp.PyArray_ZEROS(1, dims, cnp.NPY_FLOAT64, 0)
    eeg = <cnp.float64_t *>cnp.PyArray_DATA(eg)
    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

2. Test Code

The variables and functions required to reproduce the vulnerability in specfun.h and _specfun.pyx have been written into a cpp file.

#include <iostream>
#include <memory>
#include "config.h"
#include <numpy/arrayobject.h>
#include <Python.h>
enum class Status { OK = 0, NoMemory, Other};
template <typename T> Status segv(int m, int n, T c, int kd, T *cv, T *eg) {
    // =========================================================
    // Purpose: Compute the characteristic values of spheroidal
    //          wave functions
    // Input :  m  --- Mode parameter
    //          n  --- Mode parameter
    //          c  --- Spheroidal parameter
    //          KD --- Function code
    //                 KD=1 for Prolate; KD=-1 for Oblate
    // Output:  CV --- Characteristic value for given m, n and c
    //          EG(L) --- Characteristic value for mode m and n'
    //                    ( L = n' - m + 1 )
    // Return value:
    //          Status::OK
    //              Normal return.
    //          Status::NoMemory
    //              An internal memory allocation failed.
    // =========================================================

    int i, icm, j, k, k1, l, nm, nm1;
    T cs, dk0, dk1, dk2, d2k, s, t, t1, x1, xa, xb;
    // eg[<=200] is supplied by the caller

    if (c < 1e-10) {
        for (i = 1; i <= (n - m + 1); i++) {
            eg[i - 1] = (i + m) * (i + m - 1);
        }
        *cv = eg[n - m];
        return Status::OK;
    }

    // TODO: Following array sizes should be decided dynamically
    auto a = std::unique_ptr<T[]>{new (std::nothrow) T[300]()};
    auto b = std::unique_ptr<T[]>{new (std::nothrow) T[100]()};
    auto cv0 = std::unique_ptr<T[]>{new (std::nothrow) T[100]()};
    auto d = std::unique_ptr<T[]>{new (std::nothrow) T[300]()};
    auto e = std::unique_ptr<T[]>{new (std::nothrow) T[300]()};
    auto f = std::unique_ptr<T[]>{new (std::nothrow) T[300]()};
    auto g = std::unique_ptr<T[]>{new (std::nothrow) T[300]()};
    auto h = std::unique_ptr<T[]>{new (std::nothrow) T[100]()};

    if (a.get() == nullptr || b.get() == nullptr || cv0.get() == nullptr || d.get() == nullptr ||
        e.get() == nullptr || f.get() == nullptr || g.get() == nullptr || h.get() == nullptr) {
        return Status::NoMemory;
    }
    icm = (n - m + 2) / 2;
    nm = 10 + (int)(0.5 * (n - m) + c);
    cs = c * c * kd;
    k = 0;
    for (l = 0; l <= 1; l++) {

        for (i = 1; i <= nm; i++) {
            k = (l == 0 ? 2 * (i - 1) : 2 * i - 1);
            dk0 = m + k;
            dk1 = m + k + 1;
            dk2 = 2 * (m + k);
            d2k = 2 * m + k;
            a[i - 1] = (d2k + 2.0) * (d2k + 1.0) / ((dk2 + 3.0) * (dk2 + 5.0)) * cs;
            d[i - 1] = dk0 * dk1 + (2.0 * dk0 * dk1 - 2.0 * m * m - 1.0) / ((dk2 - 1.0) * (dk2 + 3.0)) * cs;
            g[i - 1] = k * (k - 1.0) / ((dk2 - 3.0) * (dk2 - 1.0)) * cs;
        }
        for (k = 2; k <= nm; k++) {
            e[k - 1] = sqrt(a[k - 2] * g[k - 1]);
            f[k - 1] = e[k - 1] * e[k - 1];
        }
        f[0] = 0.0;
        e[0] = 0.0;
        xa = d[nm - 1] + fabs(e[nm - 1]);
        xb = d[nm - 1] - fabs(e[nm - 1]);
        nm1 = nm - 1;
        for (i = 1; i <= nm1; i++) {
            t = fabs(e[i - 1]) + fabs(e[i]);
            t1 = d[i - 1] + t;
            if (xa < t1) {
                xa = t1;
            }
            t1 = d[i - 1] - t;
            if (t1 < xb) {
                xb = t1;
            }
        }
        for (i = 1; i <= icm; i++) {
            b[i - 1] = xa;
            h[i - 1] = xb;
        }
        for (k = 1; k <= icm; k++) {
            for (k1 = k; k1 <= icm; k1++) {
                if (b[k1 - 1] < b[k - 1]) {
                    b[k - 1] = b[k1 - 1];
                    break;
                }
            }
            if (k != 1) {
                if (h[k - 1] < h[k - 2]) {
                    h[k - 1] = h[k - 2];
                }
            }
            while (1) {
                x1 = (b[k - 1] + h[k - 1]) / 2.0;
                cv0[k - 1] = x1;
                if (fabs((b[k - 1] - h[k - 1]) / x1) < 1e-14) {
                    break;
                }
                j = 0;
                s = 1.0;
                for (i = 1; i <= nm; i++) {
                    if (s == 0.0) {
                        s += 1e-30;
                    }
                    t = f[i - 1] / s;
                    s = d[i - 1] - t - x1;
                    if (s < 0.0) {
                        j += 1;
                    }
                }
                if (j < k) {
                    h[k - 1] = x1;
                } else {
                    b[k - 1] = x1;
                    if (j >= icm) {
                        b[icm - 1] = x1;
                    } else {
                        if (h[j] < x1) {
                            h[j] = x1;
                        }
                        if (x1 < b[j - 1]) {
                            b[j - 1] = x1;
                        }
                    }
                }
            }
            cv0[k - 1] = x1;
            if (l == 0)
                eg[2 * k - 2] = cv0[k - 1];
            if (l == 1)
                eg[2 * k - 1] = cv0[k - 1];
        }
    }
    *cv = eg[n - m];
    return Status::OK;
}
int main () {
    Py_Initialize();
    import_array();
    int m = 1, n = 3, kd=-1;
    double c = 2.0;
    double cv;     
    npy_intp dims[1];
    dims[0] = n - m + 1;
    PyObject *eg = PyArray_Zeros(1, dims, PyArray_DescrFromType(NPY_DOUBLE), 0);
    double *eeg = (double *)PyArray_DATA((PyArrayObject *)eg);
    Status result = segv(m, n, c, kd, &cv, eeg);
    switch (result) {
        case Status::OK:
            std::cout << "cv = " << cv << std::endl;
            for (int i = 0; i < dims[0]; i++) {
                std::cout << "eeg[" << i << "] = " << eeg[i] << std::endl;
            }
            std::cout << "segv return OK" << std::endl;
            break;
        case Status::NoMemory:
            std::cout << "segv return NoMemory" << std::endl;
            break;
        case Status::Other:
            std::cout << "segv return Other" << std::endl;
    }
    Py_DECREF(eg);
    Py_Finalize();
    return 0;
}

After compiling with g++ and running, you can see successful calculations but also memory-related error messages.

# Remember to modify the paths.
g++ -g -I/opt/tools/Python-3.13.9 -I/opt/tools/Python-3.13.9/Include -I.venv/lib/python3.13/site-packages/numpy/_core/include segv.cpp /opt/tools/Python-3.13.9/lib/libpython3.13.so -o segv
./segv

cv = 10.1632
eeg[0] = 1.11855
eeg[1] = 4.22275
eeg[2] = 10.1632
segv return OK
corrupted size vs. prev_size
Aborted (core dumped)

Modify the code at lines 6058-6061 to add debug information when the length exceeds n - m + 1.

#include <iostream>
#include <memory>
#include "config.h"
#include <numpy/arrayobject.h>
#include <Python.h>
enum class Status { OK = 0, NoMemory, Other};
template <typename T> Status segv(int m, int n, T c, int kd, T *cv, T *eg) {
    int max_eg_size = n - m + 1; 
    // =========================================================
    // Purpose: Compute the characteristic values of spheroidal
    //          wave functions
    // Input :  m  --- Mode parameter
    //          n  --- Mode parameter
    //          c  --- Spheroidal parameter
    //          KD --- Function code
    //                 KD=1 for Prolate; KD=-1 for Oblate
    // Output:  CV --- Characteristic value for given m, n and c
    //          EG(L) --- Characteristic value for mode m and n'
    //                    ( L = n' - m + 1 )
    // Return value:
    //          Status::OK
    //              Normal return.
    //          Status::NoMemory
    //              An internal memory allocation failed.
    // =========================================================

    int i, icm, j, k, k1, l, nm, nm1;
    T cs, dk0, dk1, dk2, d2k, s, t, t1, x1, xa, xb;
    // eg[<=200] is supplied by the caller

    if (c < 1e-10) {
        for (i = 1; i <= (n - m + 1); i++) {
            eg[i - 1] = (i + m) * (i + m - 1);
        }
        *cv = eg[n - m];
        return Status::OK;
    }

    // TODO: Following array sizes should be decided dynamically
    auto a = std::unique_ptr<T[]>{new (std::nothrow) T[300]()};
    auto b = std::unique_ptr<T[]>{new (std::nothrow) T[100]()};
    auto cv0 = std::unique_ptr<T[]>{new (std::nothrow) T[100]()};
    auto d = std::unique_ptr<T[]>{new (std::nothrow) T[300]()};
    auto e = std::unique_ptr<T[]>{new (std::nothrow) T[300]()};
    auto f = std::unique_ptr<T[]>{new (std::nothrow) T[300]()};
    auto g = std::unique_ptr<T[]>{new (std::nothrow) T[300]()};
    auto h = std::unique_ptr<T[]>{new (std::nothrow) T[100]()};

    if (a.get() == nullptr || b.get() == nullptr || cv0.get() == nullptr || d.get() == nullptr ||
        e.get() == nullptr || f.get() == nullptr || g.get() == nullptr || h.get() == nullptr) {
        return Status::NoMemory;
    }
    icm = (n - m + 2) / 2;
    nm = 10 + (int)(0.5 * (n - m) + c);
    cs = c * c * kd;
    k = 0;
    for (l = 0; l <= 1; l++) {

        for (i = 1; i <= nm; i++) {
            k = (l == 0 ? 2 * (i - 1) : 2 * i - 1);
            dk0 = m + k;
            dk1 = m + k + 1;
            dk2 = 2 * (m + k);
            d2k = 2 * m + k;
            a[i - 1] = (d2k + 2.0) * (d2k + 1.0) / ((dk2 + 3.0) * (dk2 + 5.0)) * cs;
            d[i - 1] = dk0 * dk1 + (2.0 * dk0 * dk1 - 2.0 * m * m - 1.0) / ((dk2 - 1.0) * (dk2 + 3.0)) * cs;
            g[i - 1] = k * (k - 1.0) / ((dk2 - 3.0) * (dk2 - 1.0)) * cs;
        }
        for (k = 2; k <= nm; k++) {
            e[k - 1] = sqrt(a[k - 2] * g[k - 1]);
            f[k - 1] = e[k - 1] * e[k - 1];
        }
        f[0] = 0.0;
        e[0] = 0.0;
        xa = d[nm - 1] + fabs(e[nm - 1]);
        xb = d[nm - 1] - fabs(e[nm - 1]);
        nm1 = nm - 1;
        for (i = 1; i <= nm1; i++) {
            t = fabs(e[i - 1]) + fabs(e[i]);
            t1 = d[i - 1] + t;
            if (xa < t1) {
                xa = t1;
            }
            t1 = d[i - 1] - t;
            if (t1 < xb) {
                xb = t1;
            }
        }
        for (i = 1; i <= icm; i++) {
            b[i - 1] = xa;
            h[i - 1] = xb;
        }
        for (k = 1; k <= icm; k++) {
            for (k1 = k; k1 <= icm; k1++) {
                if (b[k1 - 1] < b[k - 1]) {
                    b[k - 1] = b[k1 - 1];
                    break;
                }
            }
            if (k != 1) {
                if (h[k - 1] < h[k - 2]) {
                    h[k - 1] = h[k - 2];
                }
            }
            while (1) {
                x1 = (b[k - 1] + h[k - 1]) / 2.0;
                cv0[k - 1] = x1;
                if (fabs((b[k - 1] - h[k - 1]) / x1) < 1e-14) {
                    break;
                }
                j = 0;
                s = 1.0;
                for (i = 1; i <= nm; i++) {
                    if (s == 0.0) {
                        s += 1e-30;
                    }
                    t = f[i - 1] / s;
                    s = d[i - 1] - t - x1;
                    if (s < 0.0) {
                        j += 1;
                    }
                }
                if (j < k) {
                    h[k - 1] = x1;
                } else {
                    b[k - 1] = x1;
                    if (j >= icm) {
                        b[icm - 1] = x1;
                    } else {
                        if (h[j] < x1) {
                            h[j] = x1;
                        }
                        if (x1 < b[j - 1]) {
                            b[j - 1] = x1;
                        }
                    }
                }
            }
            cv0[k - 1] = x1;
            int index = 0;
            if (l == 0)
                index = 2 * k - 2;
            if (l == 1)
                index = 2 * k - 1;
            if (index >= max_eg_size) {
                 std::cout << "ERROR: eg index=" << index << ", max=" << max_eg_size << std::endl;
                return Status::Other;
            } 
            eg[index] = cv0[k - 1];
        }
    }
    *cv = eg[n - m];
    return Status::OK;
}
int main () {
    Py_Initialize();
    import_array();
    int m = 1, n = 3, kd=-1;
    double c = 2.0;
    double cv;     
    npy_intp dims[1];
    dims[0] = n - m + 1;
    PyObject *eg = PyArray_Zeros(1, dims, PyArray_DescrFromType(NPY_DOUBLE), 0);
    double *eeg = (double *)PyArray_DATA((PyArrayObject *)eg);
    Status result = segv(m, n, c, kd, &cv, eeg);
    switch (result) {
        case Status::OK:
            std::cout << "cv = " << cv << std::endl;
            for (int i = 0; i < dims[0]; i++) {
                std::cout << "eeg[" << i << "] = " << eeg[i] << std::endl;
            }
            
            std::cout << "segv return OK" << std::endl;
            break;
        case Status::NoMemory:
            std::cout << "segv return NoMemory" << std::endl;
            break;
        case Status::Other:
            std::cout << "segv return Other" << std::endl;
    }
    Py_DECREF(eg);
    Py_Finalize();
    return 0;
}

After compiling with g++ and running, it can be observed that the array length is 3, but the index also reaches 3, exceeding the bounds by 1 unit.

# Remember to modify the paths.
g++ -g -I/opt/tools/Python-3.13.9 -I/opt/tools/Python-3.13.9/Include -I.venv/lib/python3.13/site-packages/numpy/_core/include segv.cpp /opt/tools/Python-3.13.9/lib/libpython3.13.so -o segv
./segv

ERROR: eg index=3, max=3
segv return Other