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 { 53 if math.Abs(real(x)) > 1 { 54 return complex(math.Pi/2, 0) // DOMAIN error 55 } 56 return complex(math.Asin(real(x)), 0) 57 } 58 ct := complex(-imag(x), real(x)) // i * x 59 xx := x * x 60 x1 := complex(1-real(xx), -imag(xx)) // 1 - x*x 61 x2 := Sqrt(x1) // x2 = sqrt(1 - x*x) 62 w := Log(ct + x2) 63 return complex(imag(w), -real(w)) // -i * w 64 } 65 66 // Asinh returns the inverse hyperbolic sine of x. 67 func Asinh(x complex128) complex128 { 68 // TODO check range 69 if imag(x) == 0 { 70 if math.Abs(real(x)) > 1 { 71 return complex(math.Pi/2, 0) // DOMAIN error 72 } 73 return complex(math.Asinh(real(x)), 0) 74 } 75 xx := x * x 76 x1 := complex(1+real(xx), imag(xx)) // 1 + x*x 77 return Log(x + Sqrt(x1)) // log(x + sqrt(1 + x*x)) 78 } 79 80 // Complex circular arc cosine 81 // 82 // DESCRIPTION: 83 // 84 // w = arccos z = PI/2 - arcsin z. 85 // 86 // ACCURACY: 87 // 88 // Relative error: 89 // arithmetic domain # trials peak rms 90 // DEC -10,+10 5200 1.6e-15 2.8e-16 91 // IEEE -10,+10 30000 1.8e-14 2.2e-15 92 93 // Acos returns the inverse cosine of x. 94 func Acos(x complex128) complex128 { 95 w := Asin(x) 96 return complex(math.Pi/2-real(w), -imag(w)) 97 } 98 99 // Acosh returns the inverse hyperbolic cosine of x. 100 func Acosh(x complex128) complex128 { 101 w := Acos(x) 102 if imag(w) <= 0 { 103 return complex(-imag(w), real(w)) // i * w 104 } 105 return complex(imag(w), -real(w)) // -i * w 106 } 107 108 // Complex circular arc tangent 109 // 110 // DESCRIPTION: 111 // 112 // If 113 // z = x + iy, 114 // 115 // then 116 // 1 ( 2x ) 117 // Re w = - arctan(-----------) + k PI 118 // 2 ( 2 2) 119 // (1 - x - y ) 120 // 121 // ( 2 2) 122 // 1 (x + (y+1) ) 123 // Im w = - log(------------) 124 // 4 ( 2 2) 125 // (x + (y-1) ) 126 // 127 // Where k is an arbitrary integer. 128 // 129 // catan(z) = -i catanh(iz). 130 // 131 // ACCURACY: 132 // 133 // Relative error: 134 // arithmetic domain # trials peak rms 135 // DEC -10,+10 5900 1.3e-16 7.8e-18 136 // IEEE -10,+10 30000 2.3e-15 8.5e-17 137 // The check catan( ctan(z) ) = z, with |x| and |y| < PI/2, 138 // had peak relative error 1.5e-16, rms relative error 139 // 2.9e-17. See also clog(). 140 141 // Atan returns the inverse tangent of x. 142 func Atan(x complex128) complex128 { 143 if real(x) == 0 && imag(x) > 1 { 144 return NaN() 145 } 146 147 x2 := real(x) * real(x) 148 a := 1 - x2 - imag(x)*imag(x) 149 if a == 0 { 150 return NaN() 151 } 152 t := 0.5 * math.Atan2(2*real(x), a) 153 w := reducePi(t) 154 155 t = imag(x) - 1 156 b := x2 + t*t 157 if b == 0 { 158 return NaN() 159 } 160 t = imag(x) + 1 161 c := (x2 + t*t) / b 162 return complex(w, 0.25*math.Log(c)) 163 } 164 165 // Atanh returns the inverse hyperbolic tangent of x. 166 func Atanh(x complex128) complex128 { 167 z := complex(-imag(x), real(x)) // z = i * x 168 z = Atan(z) 169 return complex(imag(z), -real(z)) // z = -i * z 170 } 171
