Numpy (faster?) version of @le_m's answer:

(I doubt that it is necessary to address zero division as Jyaif points out as the probability that the norm of g is 0 is mathmatically 0 and numerically almost 0, but I have added nan_to_num for now.)

from collections.abc import Sequence

import numpy as np
from numpy.typing import NDArray


def random_hypersphere(
    shape: Sequence[int], dim: int, *, surface: bool = False
) -> NDArray[np.float64]:
    r"""
    Generate random points on / in a unit hypersphere.

    Parameters
    ----------
    shape : Sequence[int]
        The shape of the output array.
    dim : int
        The dimension of the hypersphere.
    surface : bool, optional
        Whether to generate points on the surface of the hypersphere,
        by default False.

    Returns
    -------
    NDArray[np.float64]
        The generated points.

    References
    ----------
        Barthe, F., Guedon, O., Mendelson, S., & Naor, A. (2005).
        A probabilistic approach to the geometry of the \ell_p^n-ball.
        arXiv, math/0503650. Retrieved from https://arxiv.org/abs/math/0503650v1

    """
    g = np.random.normal(loc=0, scale=np.sqrt(1 / 2), size=(dim, *shape))
    if surface:
        result = g / np.linalg.vector_norm(g, axis=0)
    else:
        z = np.random.exponential(scale=1, size=shape)[None, ...]
        result = g / np.sqrt(np.sum(g**2, axis=0, keepdims=True) + z)
    return result.nan_to_num(nan=0.0)  # not sure if nan_to_num is necessary