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 tangent 32 // 33 // DESCRIPTION: 34 // 35 // If 36 // z = x + iy, 37 // 38 // then 39 // 40 // sin 2x + i sinh 2y 41 // w = --------------------. 42 // cos 2x + cosh 2y 43 // 44 // On the real axis the denominator is zero at odd multiples 45 // of PI/2. The denominator is evaluated by its Taylor 46 // series near these points. 47 // 48 // ctan(z) = -i ctanh(iz). 49 // 50 // ACCURACY: 51 // 52 // Relative error: 53 // arithmetic domain # trials peak rms 54 // DEC -10,+10 5200 7.1e-17 1.6e-17 55 // IEEE -10,+10 30000 7.2e-16 1.2e-16 56 // Also tested by ctan * ccot = 1 and catan(ctan(z)) = z. 57 58 // Tan returns the tangent of x. 59 func Tan(x complex128) complex128 { 60 d := math.Cos(2*real(x)) + math.Cosh(2*imag(x)) 61 if math.Abs(d) < 0.25 { 62 d = tanSeries(x) 63 } 64 if d == 0 { 65 return Inf() 66 } 67 return complex(math.Sin(2*real(x))/d, math.Sinh(2*imag(x))/d) 68 } 69 70 // Complex hyperbolic tangent 71 // 72 // DESCRIPTION: 73 // 74 // tanh z = (sinh 2x + i sin 2y) / (cosh 2x + cos 2y) . 75 // 76 // ACCURACY: 77 // 78 // Relative error: 79 // arithmetic domain # trials peak rms 80 // IEEE -10,+10 30000 1.7e-14 2.4e-16 81 82 // Tanh returns the hyperbolic tangent of x. 83 func Tanh(x complex128) complex128 { 84 d := math.Cosh(2*real(x)) + math.Cos(2*imag(x)) 85 if d == 0 { 86 return Inf() 87 } 88 return complex(math.Sinh(2*real(x))/d, math.Sin(2*imag(x))/d) 89 } 90 91 // Program to subtract nearest integer multiple of PI 92 func reducePi(x float64) float64 { 93 const ( 94 // extended precision value of PI: 95 DP1 = 3.14159265160560607910E0 // ?? 0x400921fb54000000 96 DP2 = 1.98418714791870343106E-9 // ?? 0x3e210b4610000000 97 DP3 = 1.14423774522196636802E-17 // ?? 0x3c6a62633145c06e 98 ) 99 t := x / math.Pi 100 if t >= 0 { 101 t += 0.5 102 } else { 103 t -= 0.5 104 } 105 t = float64(int64(t)) // int64(t) = the multiple 106 return ((x - t*DP1) - t*DP2) - t*DP3 107 } 108 109 // Taylor series expansion for cosh(2y) - cos(2x) 110 func tanSeries(z complex128) float64 { 111 const MACHEP = 1.0 / (1 << 53) 112 x := math.Abs(2 * real(z)) 113 y := math.Abs(2 * imag(z)) 114 x = reducePi(x) 115 x = x * x 116 y = y * y 117 x2 := 1.0 118 y2 := 1.0 119 f := 1.0 120 rn := 0.0 121 d := 0.0 122 for { 123 rn++ 124 f *= rn 125 rn++ 126 f *= rn 127 x2 *= x 128 y2 *= y 129 t := y2 + x2 130 t /= f 131 d += t 132 133 rn++ 134 f *= rn 135 rn++ 136 f *= rn 137 x2 *= x 138 y2 *= y 139 t = y2 - x2 140 t /= f 141 d += t 142 if !(math.Abs(t/d) > MACHEP) { 143 // Caution: Use ! and > instead of <= for correct behavior if t/d is NaN. 144 // See issue 17577. 145 break 146 } 147 } 148 return d 149 } 150 151 // Complex circular cotangent 152 // 153 // DESCRIPTION: 154 // 155 // If 156 // z = x + iy, 157 // 158 // then 159 // 160 // sin 2x - i sinh 2y 161 // w = --------------------. 162 // cosh 2y - cos 2x 163 // 164 // On the real axis, the denominator has zeros at even 165 // multiples of PI/2. Near these points it is evaluated 166 // by a Taylor series. 167 // 168 // ACCURACY: 169 // 170 // Relative error: 171 // arithmetic domain # trials peak rms 172 // DEC -10,+10 3000 6.5e-17 1.6e-17 173 // IEEE -10,+10 30000 9.2e-16 1.2e-16 174 // Also tested by ctan * ccot = 1 + i0. 175 176 // Cot returns the cotangent of x. 177 func Cot(x complex128) complex128 { 178 d := math.Cosh(2*imag(x)) - math.Cos(2*real(x)) 179 if math.Abs(d) < 0.25 { 180 d = tanSeries(x) 181 } 182 if d == 0 { 183 return Inf() 184 } 185 return complex(math.Sin(2*real(x))/d, -math.Sinh(2*imag(x))/d) 186 } 187
