diff --git a/scipy-stubs/_typing.pyi b/scipy-stubs/_typing.pyi
index da5fc9f3..0fdf1bd2 100644
--- a/scipy-stubs/_typing.pyi
+++ b/scipy-stubs/_typing.pyi
@@ -1,12 +1,22 @@
 # Helper types for internal use (type-check only).
-
 from typing import TypeAlias
 
 import numpy as np
+import numpy.typing as npt
+from typing_extensions import TypeVar
 
-__all__ = ("ArrayLike0D", "Untyped")
+__all__ = "AnyBool", "AnyChar", "AnyComplex", "AnyInt", "AnyReal", "AnyScalar", "Array0D", "Untyped"
 
+# placeholder for missing annotations
 Untyped: TypeAlias = object
 
-# keep in sync with `numpy.typing._scalar_like._ScalarLike`
-ArrayLike0D: TypeAlias = bool | int | float | complex | str | bytes | np.generic
+_SCT = TypeVar("_SCT", bound=np.generic, default=np.generic)
+Array0D: TypeAlias = np.ndarray[tuple[()], np.dtype[_SCT]]
+
+# keep in sync with `numpy._typing._scalars`
+AnyBool: TypeAlias = bool | np.bool_
+AnyInt: TypeAlias = int | np.integer[npt.NBitBase] | np.bool_
+AnyReal: TypeAlias = int | float | np.floating[npt.NBitBase] | np.integer[npt.NBitBase] | np.bool_
+AnyComplex: TypeAlias = int | float | complex | np.number[npt.NBitBase] | np.bool_
+AnyChar: TypeAlias = str | bytes
+AnyScalar: TypeAlias = int | float | complex | AnyChar | np.generic
diff --git a/scipy-stubs/spatial/_ckdtree.pyi b/scipy-stubs/spatial/_ckdtree.pyi
index d1df0aef..5f974df2 100644
--- a/scipy-stubs/spatial/_ckdtree.pyi
+++ b/scipy-stubs/spatial/_ckdtree.pyi
@@ -3,7 +3,7 @@ from typing import Any, Generic, Literal, TypeAlias, TypeVar, overload
 import numpy as np
 import numpy.typing as npt
 
-from scipy._typing import ArrayLike0D
+from scipy._typing import AnyScalar
 from scipy.sparse import coo_matrix, dok_matrix
 
 __all__ = ["cKDTree"]
@@ -123,7 +123,7 @@ class cKDTree(Generic[_BoxT]):
     def count_neighbors(
         self,
         other: cKDTree[npt.NDArray[np.float64] | None],
-        r: ArrayLike0D,
+        r: AnyScalar,
         p: float = ...,
         weights: None | tuple[None, None] = ...,
         cumulative: bool = ...,
@@ -132,7 +132,7 @@ class cKDTree(Generic[_BoxT]):
     def count_neighbors(
         self,
         other: cKDTree[npt.NDArray[np.float64] | None],
-        r: ArrayLike0D,
+        r: AnyScalar,
         p: float = ...,
         weights: _Weights = ...,
         cumulative: bool = ...,
diff --git a/scipy-stubs/special/__init__.pyi b/scipy-stubs/special/__init__.pyi
index 9d32d676..f7580e7f 100644
--- a/scipy-stubs/special/__init__.pyi
+++ b/scipy-stubs/special/__init__.pyi
@@ -1,395 +1,40 @@
+from . import _basic, _ellip_harm, _lambertw, _logsumexp, _orthogonal, _sf_error, _spfun_stats, _spherical_bessel, _ufuncs
 from ._basic import *
-from ._ellip_harm import ellip_harm, ellip_harm_2, ellip_normal
-from ._lambertw import lambertw
-from ._logsumexp import log_softmax, logsumexp, softmax
-from ._multiufuncs import *
+from ._ellip_harm import *
+from ._lambertw import *
+from ._logsumexp import *
 from ._orthogonal import *
-from ._sf_error import SpecialFunctionError, SpecialFunctionWarning
-from ._spfun_stats import multigammaln
-from ._spherical_bessel import spherical_in, spherical_jn, spherical_kn, spherical_yn
-from ._support_alternative_backends import (
-    betainc,
-    betaincc,
-    chdtr,
-    chdtrc,
-    entr,
-    erf,
-    erfc,
-    expit,
-    gammainc,
-    gammaincc,
-    gammaln,
-    i0,
-    i0e,
-    i1,
-    i1e,
-    log_ndtr,
-    logit,
-    ndtr,
-    ndtri,
-    rel_entr,
-    stdtr,
-    xlogy,
-)
+from ._sf_error import *
+from ._spfun_stats import *
+from ._spherical_bessel import *
 from ._ufuncs import *
 
-__all__ = [
-    "SpecialFunctionError",
-    "SpecialFunctionWarning",
-    "agm",
-    "ai_zeros",
-    "airy",
-    "airye",
-    "assoc_laguerre",
-    "bdtr",
-    "bdtrc",
-    "bdtri",
-    "bdtrik",
-    "bdtrin",
-    "bei",
-    "bei_zeros",
-    "beip",
-    "beip_zeros",
-    "ber",
-    "ber_zeros",
-    "bernoulli",
-    "berp",
-    "berp_zeros",
-    "besselpoly",
-    "beta",
-    "betainc",
-    "betaincc",
-    "betainccinv",
-    "betaincinv",
-    "betaln",
-    "bi_zeros",
-    "binom",
-    "boxcox",
-    "boxcox1p",
-    "btdtr",
-    "btdtri",
-    "btdtria",
-    "btdtrib",
+__all__: list[str] = [
     "c_roots",
-    "cbrt",
     "cg_roots",
-    "chdtr",
-    "chdtrc",
-    "chdtri",
-    "chdtriv",
-    "chebyc",
-    "chebys",
-    "chebyt",
-    "chebyu",
-    "chndtr",
-    "chndtridf",
-    "chndtrinc",
-    "chndtrix",
-    "clpmn",
-    "comb",
-    "cosdg",
-    "cosm1",
-    "cotdg",
-    "dawsn",
-    "digamma",
-    "diric",
-    "ellip_harm",
-    "ellip_harm_2",
-    "ellip_normal",
-    "ellipe",
-    "ellipeinc",
-    "ellipj",
-    "ellipk",
-    "ellipkinc",
-    "ellipkm1",
-    "elliprc",
-    "elliprd",
-    "elliprf",
-    "elliprg",
-    "elliprj",
-    "entr",
-    "erf",
-    "erf_zeros",
-    "erfc",
-    "erfcinv",
-    "erfcx",
-    "erfi",
-    "erfinv",
-    "errstate",
-    "euler",
-    "eval_chebyc",
-    "eval_chebys",
-    "eval_chebyt",
-    "eval_chebyu",
-    "eval_gegenbauer",
-    "eval_genlaguerre",
-    "eval_hermite",
-    "eval_hermitenorm",
-    "eval_jacobi",
-    "eval_laguerre",
-    "eval_legendre",
-    "eval_sh_chebyt",
-    "eval_sh_chebyu",
-    "eval_sh_jacobi",
-    "eval_sh_legendre",
-    "exp1",
-    "exp2",
-    "exp10",
-    "expi",
-    "expit",
-    "expm1",
-    "expn",
-    "exprel",
-    "factorial",
-    "factorial2",
-    "factorialk",
-    "fdtr",
-    "fdtrc",
-    "fdtri",
-    "fdtridfd",
-    "fresnel",
-    "fresnel_zeros",
-    "fresnelc_zeros",
-    "fresnels_zeros",
-    "gamma",
-    "gammainc",
-    "gammaincc",
-    "gammainccinv",
-    "gammaincinv",
-    "gammaln",
-    "gammasgn",
-    "gdtr",
-    "gdtrc",
-    "gdtria",
-    "gdtrib",
-    "gdtrix",
-    "gegenbauer",
-    "genlaguerre",
-    "geterr",
-    "h1vp",
-    "h2vp",
     "h_roots",
-    "hankel1",
-    "hankel1e",
-    "hankel2",
-    "hankel2e",
     "he_roots",
-    "hermite",
-    "hermitenorm",
-    "huber",
-    "hyp0f1",
-    "hyp1f1",
-    "hyp2f1",
-    "hyperu",
-    "i0",
-    "i0e",
-    "i1",
-    "i1e",
-    "inv_boxcox",
-    "inv_boxcox1p",
-    "it2i0k0",
-    "it2j0y0",
-    "it2struve0",
-    "itairy",
-    "iti0k0",
-    "itj0y0",
-    "itmodstruve0",
-    "itstruve0",
-    "iv",
-    "ive",
-    "ivp",
-    "j0",
-    "j1",
     "j_roots",
-    "jacobi",
-    "jn",
-    "jn_zeros",
-    "jnjnp_zeros",
-    "jnp_zeros",
-    "jnyn_zeros",
     "js_roots",
-    "jv",
-    "jve",
-    "jvp",
-    "k0",
-    "k0e",
-    "k1",
-    "k1e",
-    "kei",
-    "kei_zeros",
-    "keip",
-    "keip_zeros",
-    "kelvin",
-    "kelvin_zeros",
-    "ker",
-    "ker_zeros",
-    "kerp",
-    "kerp_zeros",
-    "kl_div",
-    "kn",
-    "kolmogi",
-    "kolmogorov",
-    "kv",
-    "kve",
-    "kvp",
     "l_roots",
     "la_roots",
-    "laguerre",
-    "lambertw",
-    "legendre",
-    "lmbda",
-    "log1p",
-    "log_expit",
-    "log_ndtr",
-    "log_softmax",
-    "log_wright_bessel",
-    "loggamma",
-    "logit",
-    "logsumexp",
-    "lpmn",
-    "lpmv",
-    "lpn",
-    "lqmn",
-    "lqn",
-    "mathieu_a",
-    "mathieu_b",
-    "mathieu_cem",
-    "mathieu_even_coef",
-    "mathieu_modcem1",
-    "mathieu_modcem2",
-    "mathieu_modsem1",
-    "mathieu_modsem2",
-    "mathieu_odd_coef",
-    "mathieu_sem",
-    "modfresnelm",
-    "modfresnelp",
-    "modstruve",
-    "multigammaln",
-    "nbdtr",
-    "nbdtrc",
-    "nbdtri",
-    "nbdtrik",
-    "nbdtrin",
-    "ncfdtr",
-    "ncfdtri",
-    "ncfdtridfd",
-    "ncfdtridfn",
-    "ncfdtrinc",
-    "nctdtr",
-    "nctdtridf",
-    "nctdtrinc",
-    "nctdtrit",
-    "ndtr",
-    "ndtri",
-    "ndtri_exp",
-    "nrdtrimn",
-    "nrdtrisd",
-    "obl_ang1",
-    "obl_ang1_cv",
-    "obl_cv",
-    "obl_cv_seq",
-    "obl_rad1",
-    "obl_rad1_cv",
-    "obl_rad2",
-    "obl_rad2_cv",
-    "owens_t",
     "p_roots",
-    "pbdn_seq",
-    "pbdv",
-    "pbdv_seq",
-    "pbvv",
-    "pbvv_seq",
-    "pbwa",
-    "pdtr",
-    "pdtrc",
-    "pdtri",
-    "pdtrik",
-    "perm",
-    "poch",
-    "polygamma",
-    "powm1",
-    "pro_ang1",
-    "pro_ang1_cv",
-    "pro_cv",
-    "pro_cv_seq",
-    "pro_rad1",
-    "pro_rad1_cv",
-    "pro_rad2",
-    "pro_rad2_cv",
     "ps_roots",
-    "pseudo_huber",
-    "psi",
-    "radian",
-    "rel_entr",
-    "rgamma",
-    "riccati_jn",
-    "riccati_yn",
-    "roots_chebyc",
-    "roots_chebys",
-    "roots_chebyt",
-    "roots_chebyu",
-    "roots_gegenbauer",
-    "roots_genlaguerre",
-    "roots_hermite",
-    "roots_hermitenorm",
-    "roots_jacobi",
-    "roots_laguerre",
-    "roots_legendre",
-    "roots_sh_chebyt",
-    "roots_sh_chebyu",
-    "roots_sh_jacobi",
-    "roots_sh_legendre",
-    "round",
     "s_roots",
-    "seterr",
-    "sh_chebyt",
-    "sh_chebyu",
-    "sh_jacobi",
-    "sh_legendre",
-    "shichi",
-    "sici",
-    "sinc",
-    "sindg",
-    "smirnov",
-    "smirnovi",
-    "softmax",
-    "spence",
-    "sph_harm",
-    "spherical_in",
-    "spherical_jn",
-    "spherical_kn",
-    "spherical_yn",
-    "stdtr",
-    "stdtridf",
-    "stdtrit",
-    "stirling2",
-    "struve",
     "t_roots",
-    "tandg",
-    "tklmbda",
     "ts_roots",
     "u_roots",
     "us_roots",
-    "voigt_profile",
-    "wofz",
-    "wright_bessel",
-    "wrightomega",
-    "xlog1py",
-    "xlogy",
-    "y0",
-    "y0_zeros",
-    "y1",
-    "y1_zeros",
-    "y1p_zeros",
-    "yn",
-    "yn_zeros",
-    "ynp_zeros",
-    "yv",
-    "yve",
-    "yvp",
-    "zeta",
-    "zetac",
 ]
+__all__ += _basic.__all__
+__all__ += _ellip_harm.__all__
+__all__ += _lambertw.__all__
+__all__ += _logsumexp.__all__
+__all__ += _orthogonal.__all__
+__all__ += _sf_error.__all__
+__all__ += _spfun_stats.__all__
+__all__ += _spherical_bessel.__all__
+__all__ += _ufuncs.__all__
 
 p_roots = roots_legendre
 t_roots = roots_chebyt
diff --git a/scipy-stubs/special/_add_newdocs.pyi b/scipy-stubs/special/_add_newdocs.pyi
index 90117949..7c0824fa 100644
--- a/scipy-stubs/special/_add_newdocs.pyi
+++ b/scipy-stubs/special/_add_newdocs.pyi
@@ -1,6 +1,10 @@
-from scipy._typing import Untyped
+from typing import Final
 
-docdict: dict[str, str]
+from typing_extensions import LiteralString
 
-def get(name) -> Untyped: ...
-def add_newdoc(name, doc): ...
+__all__ = ()
+
+docdict: Final[dict[str, str]]
+
+def get(name: LiteralString) -> str: ...
+def add_newdoc(name: LiteralString, doc: str) -> None: ...
diff --git a/scipy-stubs/special/_ellip_harm.pyi b/scipy-stubs/special/_ellip_harm.pyi
index 94585efd..3aafe563 100644
--- a/scipy-stubs/special/_ellip_harm.pyi
+++ b/scipy-stubs/special/_ellip_harm.pyi
@@ -1,15 +1,19 @@
 from typing import Literal
 
+import numpy as np
+
+import scipy._typing as spt
+
 __all__ = ["ellip_harm", "ellip_harm_2", "ellip_normal"]
 
 def ellip_harm(
-    h2: float,
-    k2: float,
-    n: int,
-    p: float,
-    s: float,
+    h2: spt.AnyReal,
+    k2: spt.AnyReal,
+    n: spt.AnyInt,
+    p: spt.AnyReal,
+    s: spt.AnyReal,
     signm: Literal[-1, 1] = ...,
     signn: Literal[-1, 1] = ...,
-) -> float: ...
-def ellip_harm_2(h2: float, k2: float, n: int, p: int, s: float) -> float: ...
-def ellip_normal(h2: float, k2: float, n: int, p: int) -> float: ...
+) -> np.float64: ...
+def ellip_harm_2(h2: spt.AnyReal, k2: spt.AnyReal, n: spt.AnyInt, p: spt.AnyInt, s: spt.AnyReal) -> spt.Array0D[np.float64]: ...
+def ellip_normal(h2: spt.AnyReal, k2: spt.AnyReal, n: spt.AnyReal, p: spt.AnyReal) -> spt.Array0D[np.float64]: ...
diff --git a/scipy-stubs/special/_mptestutils.pyi b/scipy-stubs/special/_mptestutils.pyi
deleted file mode 100644
index 808f8c2c..00000000
--- a/scipy-stubs/special/_mptestutils.pyi
+++ /dev/null
@@ -1,72 +0,0 @@
-from scipy._typing import Untyped
-from scipy.special._testutils import assert_func_equal as assert_func_equal
-
-class Arg:
-    def __init__(self, a=..., b=..., inclusive_a: bool = True, inclusive_b: bool = True): ...
-    def values(self, n) -> Untyped: ...
-
-class FixedArg:
-    def __init__(self, values) -> None: ...
-    def values(self, n) -> Untyped: ...
-
-class ComplexArg:
-    real: Untyped
-    imag: Untyped
-    def __init__(self, a=..., b=...) -> None: ...
-    def values(self, n) -> Untyped: ...
-
-class IntArg:
-    a: Untyped
-    b: Untyped
-    def __init__(self, a: int = -1000, b: int = 1000): ...
-    def values(self, n) -> Untyped: ...
-
-def get_args(argspec, n) -> Untyped: ...
-
-class MpmathData:
-    scipy_func: Untyped
-    mpmath_func: Untyped
-    arg_spec: Untyped
-    dps: Untyped
-    prec: Untyped
-    n: Untyped
-    rtol: Untyped
-    atol: Untyped
-    ignore_inf_sign: Untyped
-    nan_ok: Untyped
-    is_complex: Untyped
-    distinguish_nan_and_inf: Untyped
-    name: Untyped
-    param_filter: Untyped
-    def __init__(
-        self,
-        scipy_func,
-        mpmath_func,
-        arg_spec,
-        name: Untyped | None = None,
-        dps: Untyped | None = None,
-        prec: Untyped | None = None,
-        n: Untyped | None = None,
-        rtol: float = 1e-07,
-        atol: float = 1e-300,
-        ignore_inf_sign: bool = False,
-        distinguish_nan_and_inf: bool = True,
-        nan_ok: bool = True,
-        param_filter: Untyped | None = None,
-    ): ...
-    def check(self) -> Untyped: ...
-
-def assert_mpmath_equal(*a, **kw): ...
-def nonfunctional_tooslow(func) -> Untyped: ...
-def mpf2float(x) -> Untyped: ...
-def mpc2complex(x) -> Untyped: ...
-def trace_args(func) -> Untyped: ...
-
-POSIX: Untyped
-
-class TimeoutError(Exception): ...
-
-def time_limited(timeout: float = 0.5, return_val=..., use_sigalrm: bool = True) -> Untyped: ...
-def exception_to_nan(func) -> Untyped: ...
-def inf_to_nan(func) -> Untyped: ...
-def mp_assert_allclose(res, std, atol: int = 0, rtol: float = 1e-17): ...
diff --git a/scipy-stubs/special/_multiufuncs.pyi b/scipy-stubs/special/_multiufuncs.pyi
deleted file mode 100644
index 225f591f..00000000
--- a/scipy-stubs/special/_multiufuncs.pyi
+++ /dev/null
@@ -1,21 +0,0 @@
-from ._gufuncs import (
-    assoc_legendre_p_all as assoc_legendre_p_all,
-    legendre_p_all as legendre_p_all,
-    sph_harm_y_all as sph_harm_y_all,
-    sph_legendre_p_all as sph_legendre_p_all,
-)
-from ._special_ufuncs import (
-    assoc_legendre_p as assoc_legendre_p,
-    legendre_p as legendre_p,
-    sph_harm_y as sph_harm_y,
-    sph_legendre_p as sph_legendre_p,
-)
-from scipy._typing import Untyped
-
-class MultiUFunc:
-    def __init__(self, ufunc_or_ufuncs, doc: Untyped | None = None, *, force_complex_output: bool = False, **default_kwargs): ...
-    @property
-    def __doc__(self) -> Untyped: ...
-    def __call__(self, *args, **kwargs) -> Untyped: ...
-
-def _(diff_n) -> Untyped: ...
diff --git a/scipy-stubs/special/_orthogonal.pyi b/scipy-stubs/special/_orthogonal.pyi
index 7a44b595..5cc752fa 100644
--- a/scipy-stubs/special/_orthogonal.pyi
+++ b/scipy-stubs/special/_orthogonal.pyi
@@ -1,13 +1,11 @@
 from collections.abc import Callable
-from typing import Any, Literal, TypeAlias, overload
+from typing import Literal, TypeAlias, overload
 
 import numpy as np
 import numpy.typing as npt
 
-from scipy._typing import ArrayLike0D
+import scipy._typing as spt
 
-_AnyInt: TypeAlias = int | np.integer[Any]
-_AnyFloat: TypeAlias = float | np.floating[Any]
 _PointsWeights: TypeAlias = tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]
 _PointsWeightsMu: TypeAlias = tuple[npt.NDArray[np.float64], npt.NDArray[np.float64], np.float64]
 
@@ -45,97 +43,69 @@ __all__ = [
 ]
 
 @overload
-def roots_jacobi(n: _AnyInt, alpha: _AnyFloat, beta: _AnyFloat) -> _PointsWeights: ...
+def roots_jacobi(n: spt.AnyInt, alpha: spt.AnyReal, beta: spt.AnyReal, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_jacobi(n: _AnyInt, alpha: _AnyFloat, beta: _AnyFloat, mu: Literal[False]) -> _PointsWeights: ...
+def roots_jacobi(n: spt.AnyInt, alpha: spt.AnyReal, beta: spt.AnyReal, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_jacobi(n: _AnyInt, alpha: _AnyFloat, beta: _AnyFloat, mu: Literal[True]) -> _PointsWeightsMu: ...
+def roots_sh_jacobi(n: spt.AnyInt, p1: spt.AnyReal, q1: spt.AnyReal, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_sh_jacobi(n: _AnyInt, p1: _AnyFloat, q1: _AnyFloat) -> _PointsWeights: ...
+def roots_sh_jacobi(n: spt.AnyInt, p1: spt.AnyReal, q1: spt.AnyReal, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_sh_jacobi(n: _AnyInt, p1: _AnyFloat, q1: _AnyFloat, mu: Literal[False]) -> _PointsWeights: ...
+def roots_genlaguerre(n: spt.AnyInt, alpha: spt.AnyReal, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_sh_jacobi(n: _AnyInt, p1: _AnyFloat, q1: _AnyFloat, mu: Literal[True]) -> _PointsWeightsMu: ...
+def roots_genlaguerre(n: spt.AnyInt, alpha: spt.AnyReal, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_genlaguerre(n: _AnyInt, alpha: _AnyFloat) -> _PointsWeights: ...
+def roots_laguerre(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_genlaguerre(n: _AnyInt, alpha: _AnyFloat, mu: Literal[False]) -> _PointsWeights: ...
+def roots_laguerre(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_genlaguerre(n: _AnyInt, alpha: _AnyFloat, mu: Literal[True]) -> _PointsWeightsMu: ...
+def roots_hermite(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_laguerre(n: _AnyInt) -> _PointsWeights: ...
+def roots_hermite(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_laguerre(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
+def roots_hermitenorm(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_laguerre(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
+def roots_hermitenorm(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_hermite(n: _AnyInt) -> _PointsWeights: ...
+def roots_gegenbauer(n: spt.AnyInt, alpha: spt.AnyReal, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_hermite(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
+def roots_gegenbauer(n: spt.AnyInt, alpha: spt.AnyReal, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_hermite(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
+def roots_chebyt(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_hermitenorm(n: _AnyInt) -> _PointsWeights: ...
+def roots_chebyt(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_hermitenorm(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
+def roots_chebyu(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_hermitenorm(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
+def roots_chebyu(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_gegenbauer(n: _AnyInt, alpha: _AnyFloat) -> _PointsWeights: ...
+def roots_chebyc(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_gegenbauer(n: _AnyInt, alpha: _AnyFloat, mu: Literal[False]) -> _PointsWeights: ...
+def roots_chebyc(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_gegenbauer(n: _AnyInt, alpha: _AnyFloat, mu: Literal[True]) -> _PointsWeightsMu: ...
+def roots_chebys(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_chebyt(n: _AnyInt) -> _PointsWeights: ...
+def roots_chebys(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_chebyt(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
+def roots_sh_chebyt(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_chebyt(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
+def roots_sh_chebyt(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_chebyu(n: _AnyInt) -> _PointsWeights: ...
+def roots_sh_chebyu(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_chebyu(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
+def roots_sh_chebyu(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_chebyu(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
+def roots_legendre(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_chebyc(n: _AnyInt) -> _PointsWeights: ...
+def roots_legendre(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 @overload
-def roots_chebyc(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
+def roots_sh_legendre(n: spt.AnyInt, mu: Literal[False] = ...) -> _PointsWeights: ...
 @overload
-def roots_chebyc(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
-@overload
-def roots_chebys(n: _AnyInt) -> _PointsWeights: ...
-@overload
-def roots_chebys(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
-@overload
-def roots_chebys(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
-@overload
-def roots_sh_chebyt(n: _AnyInt) -> _PointsWeights: ...
-@overload
-def roots_sh_chebyt(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
-@overload
-def roots_sh_chebyt(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
-@overload
-def roots_sh_chebyu(n: _AnyInt) -> _PointsWeights: ...
-@overload
-def roots_sh_chebyu(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
-@overload
-def roots_sh_chebyu(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
-@overload
-def roots_legendre(n: _AnyInt) -> _PointsWeights: ...
-@overload
-def roots_legendre(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
-@overload
-def roots_legendre(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
-@overload
-def roots_sh_legendre(n: _AnyInt) -> _PointsWeights: ...
-@overload
-def roots_sh_legendre(n: _AnyInt, mu: Literal[False]) -> _PointsWeights: ...
-@overload
-def roots_sh_legendre(n: _AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
+def roots_sh_legendre(n: spt.AnyInt, mu: Literal[True]) -> _PointsWeightsMu: ...
 
 class orthopoly1d(np.poly1d):
+    @property
+    def limits(self) -> tuple[float, float]: ...
     def __init__(
         self,
         roots: npt.ArrayLike,
@@ -147,32 +117,24 @@ class orthopoly1d(np.poly1d):
         monic: bool = ...,
         eval_func: np.ufunc = ...,
     ) -> None: ...
-    @property
-    def limits(self) -> tuple[float, float]: ...
     def weight_func(self, x: float) -> float: ...
-    @overload
-    def __call__(self, x: ArrayLike0D) -> Any: ...
-    @overload
-    def __call__(self, x: np.poly1d) -> np.poly1d: ...
-    @overload
-    def __call__(self, x: npt.ArrayLike) -> npt.NDArray[np.float64] | npt.NDArray[np.complex128] | npt.NDArray[np.object_]: ...
 
-def legendre(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def chebyt(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def chebyu(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def chebyc(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def chebys(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def jacobi(n: _AnyInt, alpha: _AnyFloat, beta: _AnyFloat, monic: bool = ...) -> orthopoly1d: ...
-def laguerre(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def genlaguerre(n: _AnyInt, alpha: _AnyFloat, monic: bool = ...) -> orthopoly1d: ...
-def hermite(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def hermitenorm(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def gegenbauer(n: _AnyInt, alpha: _AnyFloat, monic: bool = ...) -> orthopoly1d: ...
-def sh_legendre(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def sh_chebyt(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def sh_chebyu(n: _AnyInt, monic: bool = ...) -> orthopoly1d: ...
-def sh_jacobi(n: _AnyInt, p: _AnyFloat, q: _AnyFloat, monic: bool = ...) -> orthopoly1d: ...
+def legendre(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def chebyt(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def chebyu(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def chebyc(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def chebys(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def jacobi(n: spt.AnyInt, alpha: spt.AnyReal, beta: spt.AnyReal, monic: bool = ...) -> orthopoly1d: ...
+def laguerre(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def genlaguerre(n: spt.AnyInt, alpha: spt.AnyReal, monic: bool = ...) -> orthopoly1d: ...
+def hermite(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def hermitenorm(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def gegenbauer(n: spt.AnyInt, alpha: spt.AnyReal, monic: bool = ...) -> orthopoly1d: ...
+def sh_legendre(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def sh_chebyt(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def sh_chebyu(n: spt.AnyInt, monic: bool = ...) -> orthopoly1d: ...
+def sh_jacobi(n: spt.AnyInt, p: spt.AnyReal, q: spt.AnyReal, monic: bool = ...) -> orthopoly1d: ...
 
 # These functions are not public, but still need stubs because they
 # get checked in the tests.
-def _roots_hermite_asy(n: _AnyInt) -> _PointsWeights: ...
+def _roots_hermite_asy(n: spt.AnyInt) -> _PointsWeights: ...
diff --git a/scipy-stubs/special/_test_internal.pyi b/scipy-stubs/special/_test_internal.pyi
deleted file mode 100644
index db0dfdb7..00000000
--- a/scipy-stubs/special/_test_internal.pyi
+++ /dev/null
@@ -1,9 +0,0 @@
-import numpy as np
-
-def have_fenv() -> bool: ...
-def random_double(size: int) -> np.float64: ...
-def test_add_round(size: int, mode: str) -> None: ...
-
-def _dd_exp(xhi: float, xlo: float) -> tuple[float, float]: ...
-def _dd_log(xhi: float, xlo: float) -> tuple[float, float]: ...
-def _dd_expm1(xhi: float, xlo: float) -> tuple[float, float]: ...
diff --git a/scipy-stubs/special/cython_special.pyi b/scipy-stubs/special/cython_special.pyi
index 6907d9ce..e63d5803 100644
--- a/scipy-stubs/special/cython_special.pyi
+++ b/scipy-stubs/special/cython_special.pyi
@@ -211,7 +211,7 @@ __all__ = [
 
 __pyx_capi__: Final[Mapping[str, CapsuleType]]
 
-_CythonFunctionOrMethod: TypeAlias = Callable[..., object]
+_CythonFunctionOrMethod: TypeAlias = Callable[..., object]  # no-any-explicit
 
 agm: _CythonFunctionOrMethod
 bdtr: _CythonFunctionOrMethod
diff --git a/scipy-stubs/special/utils/convert.pyi b/scipy-stubs/special/utils/convert.pyi
deleted file mode 100644
index eb6d21ad..00000000
--- a/scipy-stubs/special/utils/convert.pyi
+++ /dev/null
@@ -1,13 +0,0 @@
-from scipy._typing import Untyped
-
-DATA_DIR: str
-BOOST_SRC: str
-CXX_COMMENT: Untyped
-DATA_REGEX: Untyped
-ITEM_REGEX: Untyped
-HEADER_REGEX: Untyped
-IGNORE_PATTERNS: Untyped
-
-def parse_ipp_file(filename) -> Untyped: ...
-def dump_dataset(filename, data): ...
-def dump_datasets(filename): ...
diff --git a/scipy-stubs/special/utils/datafunc.pyi b/scipy-stubs/special/utils/datafunc.pyi
deleted file mode 100644
index 404cfb38..00000000
--- a/scipy-stubs/special/utils/datafunc.pyi
+++ /dev/null
@@ -1,4 +0,0 @@
-from scipy._typing import Untyped
-
-def parse_txt_data(filename) -> Untyped: ...
-def run_test(filename, funcs, args=[0]) -> Untyped: ...
diff --git a/scipy-stubs/special/utils/makenpz.pyi b/scipy-stubs/special/utils/makenpz.pyi
deleted file mode 100644
index fbbdad44..00000000
--- a/scipy-stubs/special/utils/makenpz.pyi
+++ /dev/null
@@ -1,4 +0,0 @@
-from scipy._typing import Untyped
-
-def newer(source, target) -> Untyped: ...
-def main(): ...