I opened gh-21404 to suggest to add a function `scipy.special.pow1p(x,y)`

to evaluate (1+x)^y accurately for x \ge -1 and y \in \mathbb{R}.

The function would be particularly useful for x close to 1, in which case the obvious formula `(1+x)**y`

loses accuracy. Such calculation is used e.g. by the probability mass function of the geometric distribution: (1-p)^k p. The function would also improve accuracy for moderate x and large y, reducing the relative error from `1e-13`

to `1e-16`

.

The same function was recently introduced in R. A similar function called `compound`

is defined in IEEE-754, with y restricted to integers.

Do you think itâ€™s worth having the `pow1p`

function in `scipy.special`

?