atanh.go 1.9 KB

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