def proj_l2(v, gamma, axis=None):
    r"""Compute the projection operator of the :math:`\ell_2` norm.

    The projection operator of the uncentered :math:`\ell_2` norm,

    .. math::
      \mathrm{argmin}_{\mathbf{x}} (1/2) \| \mathbf{x} - \mathbf{v} \|_2^2 \;
      \text{ s.t. } \; \| \mathbf{x} - \mathbf{s} \|_2 \leq \gamma

    can be computed as :math:`\mathbf{s} + \mathrm{proj}_{f,\gamma}
    (\mathbf{v} - \mathbf{s})` where :math:`f(\mathbf{x}) =
    \| \mathbf{x} \|_2`.


    Parameters
    ----------
    v : array_like
      Input array :math:`\mathbf{v}`
    gamma : float
      Parameter :math:`\gamma`
    axis : None or int or tuple of ints, optional (default None)
      Axes of `v` over which to compute the :math:`\ell_2` norm. If
      `None`, an entire multi-dimensional array is treated as a vector.
      If axes are specified, then distinct norm values are computed
      over the indices of the remaining axes of input array `v`.

    Returns
    -------
    x : ndarray
      Output array
    """

    d = np.sqrt(np.sum(v**2, axis=axis, keepdims=True))
    return np.asarray((d <= gamma) * v +
                      (d > gamma) * (gamma * sl.zdivide(v, d)),
                      dtype=v.dtype)