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