atanh.go 1.7 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  1. package math32
  2. // The original C code, the long comment, and the constants
  3. // below are from FreeBSD's /usr/src/lib/msun/src/e_atanh.c
  4. // and came with this notice. The go code is a simplified
  5. // version of the original C.
  6. //
  7. // ====================================================
  8. // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  9. //
  10. // Developed at SunPro, a Sun Microsystems, Inc. business.
  11. // Permission to use, copy, modify, and distribute this
  12. // software is freely granted, provided that this notice
  13. // is preserved.
  14. // ====================================================
  15. //
  16. //
  17. // __ieee754_atanh(x)
  18. // Method :
  19. // 1. Reduce x to positive by atanh(-x) = -atanh(x)
  20. // 2. For x>=0.5
  21. // 1 2x x
  22. // atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
  23. // 2 1 - x 1 - x
  24. //
  25. // For x<0.5
  26. // atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
  27. //
  28. // Special cases:
  29. // atanh(x) is NaN if |x| > 1 with signal;
  30. // atanh(NaN) is that NaN with no signal;
  31. // atanh(+-1) is +-INF with signal.
  32. //
  33. // Atanh returns the inverse hyperbolic tangent of x.
  34. //
  35. // Special cases are:
  36. // Atanh(1) = +Inf
  37. // Atanh(±0) = ±0
  38. // Atanh(-1) = -Inf
  39. // Atanh(x) = NaN if x < -1 or x > 1
  40. // Atanh(NaN) = NaN
  41. func Atanh(x float32) float32 {
  42. const NearZero = 1.0 / (1 << 28) // 2**-28
  43. // special cases
  44. switch {
  45. case x < -1 || x > 1 || IsNaN(x):
  46. return NaN()
  47. case x == 1:
  48. return Inf(1)
  49. case x == -1:
  50. return Inf(-1)
  51. }
  52. sign := false
  53. if x < 0 {
  54. x = -x
  55. sign = true
  56. }
  57. var temp float32
  58. switch {
  59. case x < NearZero:
  60. temp = x
  61. case x < 0.5:
  62. temp = x + x
  63. temp = 0.5 * Log1p(temp+temp*x/(1-x))
  64. default:
  65. temp = 0.5 * Log1p((x+x)/(1-x))
  66. }
  67. if sign {
  68. temp = -temp
  69. }
  70. return temp
  71. }