Optional output array for the function results

    Returns
    -------
    z : scalar or ndarray
        Computed x*log1p(y)

    Notes
    -----

    .. versionadded:: 0.13.0

    Examples
    --------
    This example shows how the function can be used to calculate the log of
    the probability mass function for a geometric discrete random variable.
    The probability mass function of the geometric distribution is defined
    as follows:

    .. math:: f(k) = (1-p)^{k-1} p

    where :math:`p` is the probability of a single success
    and :math:`1-p` is the probability of a single failure
    and :math:`k` is the number of trials to get the first success.

    >>> import numpy as np
    >>> from scipy.special import xlog1py
    >>> p = 0.5
    >>> k = 100
    >>> _pmf = np.power(1 - p, k - 1) * p
    >>> _pmf
    7.888609052210118e-31

    If we take k as a relatively large number the value of the probability
    mass function can become very low. In such cases taking the log of the
    pmf would be more suitable as the log function can change the values
    to a scale that is more appropriate to work with.

    >>> _log_pmf = xlog1py(k - 1, -p) + np.log(p)
    >>> _log_pmf
    -69.31471805599453

    We can confirm that we get a value close to the original pmf value by
    taking the exponential of the log pmf.

    >>> _orig_pmf = np.exp(_log_pmf)
    >>> np.isclose(_pmf, _orig_pmf)
    True

    