--------
    >>> import numpy as np
    >>> from scipy.special import powm1

    >>> x = np.array([1.2, 10.0, 0.9999999975])
    >>> y = np.array([1e-9, 1e-11, 0.1875])
    >>> powm1(x, y)
    array([ 1.82321557e-10,  2.30258509e-11, -4.68749998e-10])

    It can be verified that the relative errors in those results
    are less than 2.5e-16.

    Compare that to the result of ``x**y - 1``, where the
    relative errors are all larger than 8e-8:

    >>> x**y - 1
    array([ 1.82321491e-10,  2.30258035e-11, -4.68750039e-10])

    Ú