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?