atanhf.c 926 B

12345678910111213141516171819202122232425262728293031323334
  1. /*****************************************************************************/
  2. /*****************************************************************************/
  3. // atanhf from musl-0.9.15
  4. /*****************************************************************************/
  5. /*****************************************************************************/
  6. #include "libm.h"
  7. /* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
  8. float atanhf(float x)
  9. {
  10. union {float f; uint32_t i;} u = {.f = x};
  11. unsigned s = u.i >> 31;
  12. float_t y;
  13. /* |x| */
  14. u.i &= 0x7fffffff;
  15. y = u.f;
  16. if (u.i < 0x3f800000 - (1<<23)) {
  17. if (u.i < 0x3f800000 - (32<<23)) {
  18. /* handle underflow */
  19. if (u.i < (1<<23))
  20. FORCE_EVAL((float)(y*y));
  21. } else {
  22. /* |x| < 0.5, up to 1.7ulp error */
  23. y = 0.5f*log1pf(2*y + 2*y*y/(1-y));
  24. }
  25. } else {
  26. /* avoid overflow */
  27. y = 0.5f*log1pf(2*(y/(1-y)));
  28. }
  29. return s ? -y : y;
  30. }