asin.go 4.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171
  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 cmplx
  5. import "math"
  6. // The original C code, the long comment, and the constants
  7. // below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
  8. // The go code is a simplified version of the original C.
  9. //
  10. // Cephes Math Library Release 2.8: June, 2000
  11. // Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
  12. //
  13. // The readme file at http://netlib.sandia.gov/cephes/ says:
  14. // Some software in this archive may be from the book _Methods and
  15. // Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
  16. // International, 1989) or from the Cephes Mathematical Library, a
  17. // commercial product. In either event, it is copyrighted by the author.
  18. // What you see here may be used freely but it comes with no support or
  19. // guarantee.
  20. //
  21. // The two known misprints in the book are repaired here in the
  22. // source listings for the gamma function and the incomplete beta
  23. // integral.
  24. //
  25. // Stephen L. Moshier
  26. // moshier@na-net.ornl.gov
  27. // Complex circular arc sine
  28. //
  29. // DESCRIPTION:
  30. //
  31. // Inverse complex sine:
  32. // 2
  33. // w = -i clog( iz + csqrt( 1 - z ) ).
  34. //
  35. // casin(z) = -i casinh(iz)
  36. //
  37. // ACCURACY:
  38. //
  39. // Relative error:
  40. // arithmetic domain # trials peak rms
  41. // DEC -10,+10 10100 2.1e-15 3.4e-16
  42. // IEEE -10,+10 30000 2.2e-14 2.7e-15
  43. // Larger relative error can be observed for z near zero.
  44. // Also tested by csin(casin(z)) = z.
  45. // Asin returns the inverse sine of x.
  46. func Asin(x complex128) complex128 {
  47. if imag(x) == 0 {
  48. if math.Abs(real(x)) > 1 {
  49. return complex(math.Pi/2, 0) // DOMAIN error
  50. }
  51. return complex(math.Asin(real(x)), 0)
  52. }
  53. ct := complex(-imag(x), real(x)) // i * x
  54. xx := x * x
  55. x1 := complex(1-real(xx), -imag(xx)) // 1 - x*x
  56. x2 := Sqrt(x1) // x2 = sqrt(1 - x*x)
  57. w := Log(ct + x2)
  58. return complex(imag(w), -real(w)) // -i * w
  59. }
  60. // Asinh returns the inverse hyperbolic sine of x.
  61. func Asinh(x complex128) complex128 {
  62. // TODO check range
  63. if imag(x) == 0 {
  64. if math.Abs(real(x)) > 1 {
  65. return complex(math.Pi/2, 0) // DOMAIN error
  66. }
  67. return complex(math.Asinh(real(x)), 0)
  68. }
  69. xx := x * x
  70. x1 := complex(1+real(xx), imag(xx)) // 1 + x*x
  71. return Log(x + Sqrt(x1)) // log(x + sqrt(1 + x*x))
  72. }
  73. // Complex circular arc cosine
  74. //
  75. // DESCRIPTION:
  76. //
  77. // w = arccos z = PI/2 - arcsin z.
  78. //
  79. // ACCURACY:
  80. //
  81. // Relative error:
  82. // arithmetic domain # trials peak rms
  83. // DEC -10,+10 5200 1.6e-15 2.8e-16
  84. // IEEE -10,+10 30000 1.8e-14 2.2e-15
  85. // Acos returns the inverse cosine of x.
  86. func Acos(x complex128) complex128 {
  87. w := Asin(x)
  88. return complex(math.Pi/2-real(w), -imag(w))
  89. }
  90. // Acosh returns the inverse hyperbolic cosine of x.
  91. func Acosh(x complex128) complex128 {
  92. w := Acos(x)
  93. if imag(w) <= 0 {
  94. return complex(-imag(w), real(w)) // i * w
  95. }
  96. return complex(imag(w), -real(w)) // -i * w
  97. }
  98. // Complex circular arc tangent
  99. //
  100. // DESCRIPTION:
  101. //
  102. // If
  103. // z = x + iy,
  104. //
  105. // then
  106. // 1 ( 2x )
  107. // Re w = - arctan(-----------) + k PI
  108. // 2 ( 2 2)
  109. // (1 - x - y )
  110. //
  111. // ( 2 2)
  112. // 1 (x + (y+1) )
  113. // Im w = - log(------------)
  114. // 4 ( 2 2)
  115. // (x + (y-1) )
  116. //
  117. // Where k is an arbitrary integer.
  118. //
  119. // catan(z) = -i catanh(iz).
  120. //
  121. // ACCURACY:
  122. //
  123. // Relative error:
  124. // arithmetic domain # trials peak rms
  125. // DEC -10,+10 5900 1.3e-16 7.8e-18
  126. // IEEE -10,+10 30000 2.3e-15 8.5e-17
  127. // The check catan( ctan(z) ) = z, with |x| and |y| < PI/2,
  128. // had peak relative error 1.5e-16, rms relative error
  129. // 2.9e-17. See also clog().
  130. // Atan returns the inverse tangent of x.
  131. func Atan(x complex128) complex128 {
  132. if real(x) == 0 && imag(x) > 1 {
  133. return NaN()
  134. }
  135. x2 := real(x) * real(x)
  136. a := 1 - x2 - imag(x)*imag(x)
  137. if a == 0 {
  138. return NaN()
  139. }
  140. t := 0.5 * math.Atan2(2*real(x), a)
  141. w := reducePi(t)
  142. t = imag(x) - 1
  143. b := x2 + t*t
  144. if b == 0 {
  145. return NaN()
  146. }
  147. t = imag(x) + 1
  148. c := (x2 + t*t) / b
  149. return complex(w, 0.25*math.Log(c))
  150. }
  151. // Atanh returns the inverse hyperbolic tangent of x.
  152. func Atanh(x complex128) complex128 {
  153. z := complex(-imag(x), real(x)) // z = i * x
  154. z = Atan(z)
  155. return complex(imag(z), -real(z)) // z = -i * z
  156. }