Add the pow1p function

[mborland]

This implements the pow1p function as discussed and implemented in scipy/scipy#21404 and scipy/scipy#21410. Simplifies a couple of the branches and adds in policy based design for the error paths. This version is also tested with real_concept and on the GPU platforms (NVCC, NVRTC, and SYCL).

@fancidev let me know if you see any errors in simplification. The one to call out is here: https://github.com/boostorg/math/compare/pow1p?expand=1#diff-1a6b6446fd424370a1fe6e67b8166e41ab4805c1ecfc4ca15ceae7492f4023a5R137. I fall back to pow instead of throwing on non-integer power.

@NAThompson and @jzmaddock your review would be appreciated as always.

CC: @WarrenWeckesser

[fancidev]

It’s really cool to have this in boost! Thanks a ton!

There is some refinement I’m looking to apply: (1) represent the result of $s/t$ using double-double, and (2) expand $\ln(1+u)$ to second order. I think with these refinements
the implementation is (more or less) provably accurate to a few machine epsilon. The current code still has some input that leads to relative error over 1e-14, which I was able to construct by an error analysis.

I’m going to work out the proof on paper first, and then suggest the changes here.

Another optimization that could be applied is that if a floating point type is available whose mantissa has more bits than the total bits of T, such as double vs float or x87 long double vs double, then it should be accurate enough (for T) to compute $(1+x)^y$ using exp(y*log1p(x)), provided *, log1p, and exp are accurately computed using the extended precision type. I’lol try verify this on paper and through numerical experiments.

Finally, I’m going to prepare some edge input as test cases, to complement your randomized test cases, which I find really helpful!

P.S. The (floating point) exception signaling (via the policy object) might have some corner cases to refine, but I’d leave that after the computational part is finalized.

Thanks again for bring this into boost!!

[fancidev]

I made an update to the algorithm (in the referenced PR) that makes the evaluation more accurate for tiny $x$ and huge $y$:

scipy/scipy@a3f89c5

I think the method is now accurate to a few epsilon relative error for all $x > -1$ and all $y$, and I expect minimal changes to this part in the future. A bit more update may come for the case $x < -1$. I'm still in the progress of adding more test cases and completing the algorithm description.

@mborland If you could update line 193-219 of pow1p.hpp with the changes it would be great! Otherwise I'll suggest a patch off your PR after I finish the remaining points.

[mborland]

I made an update to the algorithm (in the referenced PR) that makes the evaluation more accurate for tiny x and huge y :

scipy/scipy@a3f89c5

I think the method is now accurate to a few epsilon relative error for all x > − 1 and all y , and I expect minimal changes to this part in the future. A bit more update may come for the case x < − 1 . I'm still in the progress of adding more test cases and completing the algorithm description.

@mborland If you could update line 193-219 of pow1p.hpp with the changes it would be great! Otherwise I'll suggest a patch off your PR after I finish the remaining points.

Changes have been applied.

Edit: @fancidev as you work on test cases it looks like right now I am not hitting anything in your double-T range with the random values. If you could come up with cases to exercise those it would be helpful.

[codecov]

Codecov Report

Attention: Patch coverage is 92.07921% with 8 lines in your changes missing coverage. Please review.

Project coverage is 94.06%. Comparing base (66b362a) to head (f7150ee).

Files	Patch %	Lines
include/boost/math/special_functions/pow1p.hpp	87.30%	8 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1183      +/-   ##
===========================================
- Coverage    94.08%   94.06%   -0.03%     
===========================================
  Files          780      782       +2     
  Lines        65815    65916     +101     
===========================================
+ Hits         61922    62003      +81     
- Misses        3893     3913      +20     

Files	Coverage Δ
test/test_pow1p.cpp	100.00% <100.00%> (ø)
include/boost/math/special_functions/pow1p.hpp	87.30% <87.30%> (ø)

... and 1 file with indirect coverage changes

---

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 66b362a...f7150ee. Read the comment docs.