Repository: scipy/scipy
File: scipy/special/_special_ufuncs.cpp (C++)

Code before:

    PyObject *_special_ufuncs = PyModule_Create(&_special_ufuncs_def);
    if (_special_ufuncs == nullptr) {
        return NULL;
    }

#if Py_GIL_DISABLED
    PyUnstable_Module_SetGIL(_special_ufuncs, Py_MOD_GIL_NOT_USED);
#endif

    PyObject *_cospi =
        xsf::numpy::ufunc({static_cast<xsf::numpy::f_f>(xsf::cospi), static_cast<xsf::numpy::d_d>(xsf::cospi),
                           static_cast<xsf::numpy::F_F>(xsf::cospi), static_cast<xsf::numpy::D_D>(xsf::cospi)},
                          "_cospi", _cospi_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_cospi", _cospi);

    PyObject *_lambertw = xsf::numpy::ufunc(
        {static_cast<xsf::numpy::Dld_D>(xsf::lambertw), static_cast<xsf::numpy::Flf_F>(xsf::lambertw)}, "_lambertw",
        lambertw_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_lambertw", _lambertw);

    PyObject *_scaled_exp1 = xsf::numpy::ufunc(
        {static_cast<xsf::numpy::d_d>(xsf::scaled_exp1), static_cast<xsf::numpy::f_f>(xsf::scaled_exp1)},
        "_scaled_exp1", scaled_exp1_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_scaled_exp1", _scaled_exp1);

    PyObject *_sinpi =
        xsf::numpy::ufunc({static_cast<xsf::numpy::f_f>(xsf::sinpi), static_cast<xsf::numpy::d_d>(xsf::sinpi),
                           static_cast<xsf::numpy::F_F>(xsf::sinpi), static_cast<xsf::numpy::D_D>(xsf::sinpi)},
                          "_sinpi", _sinpi_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_sinpi", _sinpi);

    PyObject *_zeta =
        xsf::numpy::ufunc({static_cast<xsf::numpy::ff_f>(xsf::zeta), static_cast<xsf::numpy::Ff_F>(xsf::zeta),
                           static_cast<xsf::numpy::dd_d>(xsf::zeta), static_cast<xsf::numpy::Dd_D>(xsf::zeta)},
                          "_zeta", _zeta_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_zeta", _zeta);

    PyObject *airy =
        xsf::numpy::ufunc({static_cast<xsf::numpy::f_ffff>(xsf::airy), static_cast<xsf::numpy::d_dddd>(xsf::airy),
                           static_cast<xsf::numpy::F_FFFF>(xsf::airy), static_cast<xsf::numpy::D_DDDD>(xsf::airy)},
                          4, "airy", airy_doc);
    PyModule_AddObjectRef(_special_ufuncs, "airy", airy);

    PyObject *airye =
        xsf::numpy::ufunc({static_cast<xsf::numpy::f_ffff>(xsf::airye), static_cast<xsf::numpy::d_dddd>(xsf::airye),
                           static_cast<xsf::numpy::F_FFFF>(xsf::airye), static_cast<xsf::numpy::D_DDDD>(xsf::airye)},
                          4, "airye", airye_doc);
    PyModule_AddObjectRef(_special_ufuncs, "airye", airye);

    PyObject *bei = xsf::numpy::ufunc({static_cast<xsf::numpy::f_f>(xsf::bei), static_cast<xsf::numpy::d_d>(xsf::bei)},

Review comment (suggestion):
Just pointing this out for the general audience. This  works because `gen_harmonic` and `normalized_gen_harmonic` each only have one type signature so we don't need to specify an overload with something like `static_cast<xsf::numpy::qd_d>(normalized_gen_harmonic)`.

 My preference would be to include the `static_cast` anyway just for the sake of consistency and readability. Technically though, `q` is for `long long` which in practice for every platform SciPy supports should be the same as `int64`, then again, who knows what the future holds. Do you think it would be better if there more explicit signatures in `xsf/numpy.h` like `xsf::numpy::i8d_d`? I'm fine leaving this as is for now.

There's also a question of whether we may actually additional overloads, but as is, this does reproduce the behavior with respect to types of existing `zipfian` functions, since the old `_gen_harmonic` explicitly made the output have `dtype` `float`. That's a question I'm not really concerned about for now, but just noting it for the future.

I'd also like to point out for the general audience that the current implementation does avoid issues surrounding NumPy default ints on different platforms (https://github.com/scipy/scipy/pull/21401) since the ufunc explicitly uses `int64`. This ties in with the question above of whether we want different overloads, since maybe it would be nice to allow computation with native 32 bit ints on 32 bit platforms, but it's not something I think matters that much for now.

Code after:
    if (PyErr_Occurred()) {
        return NULL;
    }

    PyObject *_special_ufuncs = PyModule_Create(&_special_ufuncs_def);
    if (_special_ufuncs == nullptr) {
        return NULL;
    }

#if Py_GIL_DISABLED
    PyUnstable_Module_SetGIL(_special_ufuncs, Py_MOD_GIL_NOT_USED);
#endif

    PyObject *_cospi =
        xsf::numpy::ufunc({static_cast<xsf::numpy::f_f>(xsf::cospi), static_cast<xsf::numpy::d_d>(xsf::cospi),
                           static_cast<xsf::numpy::F_F>(xsf::cospi), static_cast<xsf::numpy::D_D>(xsf::cospi)},
                          "_cospi", _cospi_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_cospi", _cospi);

    PyObject *_gen_harmonic =
        xsf::numpy::ufunc({gen_harmonic}, "_gen_harmonic", _gen_harmonic_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_gen_harmonic", _gen_harmonic);

    PyObject *_normalized_gen_harmonic =
        xsf::numpy::ufunc({normalized_gen_harmonic}, "_normalized_gen_harmonic", _normalized_gen_harmonic_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_normalized_gen_harmonic", _normalized_gen_harmonic);

    PyObject *_lambertw = xsf::numpy::ufunc(
        {static_cast<xsf::numpy::Dld_D>(xsf::lambertw), static_cast<xsf::numpy::Flf_F>(xsf::lambertw)}, "_lambertw",
        lambertw_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_lambertw", _lambertw);

    PyObject *_scaled_exp1 = xsf::numpy::ufunc(
        {static_cast<xsf::numpy::d_d>(xsf::scaled_exp1), static_cast<xsf::numpy::f_f>(xsf::scaled_exp1)},
        "_scaled_exp1", scaled_exp1_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_scaled_exp1", _scaled_exp1);

    PyObject *_sinpi =
        xsf::numpy::ufunc({static_cast<xsf::numpy::f_f>(xsf::sinpi), static_cast<xsf::numpy::d_d>(xsf::sinpi),
                           static_cast<xsf::numpy::F_F>(xsf::sinpi), static_cast<xsf::numpy::D_D>(xsf::sinpi)},
                          "_sinpi", _sinpi_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_sinpi", _sinpi);

    PyObject *_zeta =
        xsf::numpy::ufunc({static_cast<xsf::numpy::ff_f>(xsf::zeta), static_cast<xsf::numpy::Ff_F>(xsf::zeta),
                           static_cast<xsf::numpy::dd_d>(xsf::zeta), static_cast<xsf::numpy::Dd_D>(xsf::zeta)},
                          "_zeta", _zeta_doc);
    PyModule_AddObjectRef(_special_ufuncs, "_zeta", _zeta);

    PyObject *airy =
        xsf::numpy::ufunc({static_cast<xsf::numpy::f_ffff>(xsf::airy), static_cast<xsf::numpy::d_dddd>(xsf::airy),