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 ( 8 "math" 9 "math/bits" 10 ) 11 12 // The original C code, the long comment, and the constants 13 // below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c. 14 // The go code is a simplified version of the original C. 15 // 16 // Cephes Math Library Release 2.8: June, 2000 17 // Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier 18 // 19 // The readme file at http://netlib.sandia.gov/cephes/ says: 20 // Some software in this archive may be from the book _Methods and 21 // Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster 22 // International, 1989) or from the Cephes Mathematical Library, a 23 // commercial product. In either event, it is copyrighted by the author. 24 // What you see here may be used freely but it comes with no support or 25 // guarantee. 26 // 27 // The two known misprints in the book are repaired here in the 28 // source listings for the gamma function and the incomplete beta 29 // integral. 30 // 31 // Stephen L. Moshier 32 // [email protected] 33 34 // Complex circular tangent 35 // 36 // DESCRIPTION: 37 // 38 // If 39 // z = x + iy, 40 // 41 // then 42 // 43 // sin 2x + i sinh 2y 44 // w = --------------------. 45 // cos 2x + cosh 2y 46 // 47 // On the real axis the denominator is zero at odd multiples 48 // of PI/2. The denominator is evaluated by its Taylor 49 // series near these points. 50 // 51 // ctan(z) = -i ctanh(iz). 52 // 53 // ACCURACY: 54 // 55 // Relative error: 56 // arithmetic domain # trials peak rms 57 // DEC -10,+10 5200 7.1e-17 1.6e-17 58 // IEEE -10,+10 30000 7.2e-16 1.2e-16 59 // Also tested by ctan * ccot = 1 and catan(ctan(z)) = z. 60 61 // Tan returns the tangent of x. 62 func Tan(x complex128) complex128 { 63 switch re, im := real(x), imag(x); { 64 case math.IsInf(im, 0): 65 switch { 66 case math.IsInf(re, 0) || math.IsNaN(re): 67 return complex(math.Copysign(0, re), math.Copysign(1, im)) 68 } 69 return complex(math.Copysign(0, math.Sin(2*re)), math.Copysign(1, im)) 70 case re == 0 && math.IsNaN(im): 71 return x 72 } 73 d := math.Cos(2*real(x)) + math.Cosh(2*imag(x)) 74 if math.Abs(d) < 0.25 { 75 d = tanSeries(x) 76 } 77 if d == 0 { 78 return Inf() 79 } 80 return complex(math.Sin(2*real(x))/d, math.Sinh(2*imag(x))/d) 81 } 82 83 // Complex hyperbolic tangent 84 // 85 // DESCRIPTION: 86 // 87 // tanh z = (sinh 2x + i sin 2y) / (cosh 2x + cos 2y) . 88 // 89 // ACCURACY: 90 // 91 // Relative error: 92 // arithmetic domain # trials peak rms 93 // IEEE -10,+10 30000 1.7e-14 2.4e-16 94 95 // Tanh returns the hyperbolic tangent of x. 96 func Tanh(x complex128) complex128 { 97 switch re, im := real(x), imag(x); { 98 case math.IsInf(re, 0): 99 switch { 100 case math.IsInf(im, 0) || math.IsNaN(im): 101 return complex(math.Copysign(1, re), math.Copysign(0, im)) 102 } 103 return complex(math.Copysign(1, re), math.Copysign(0, math.Sin(2*im))) 104 case im == 0 && math.IsNaN(re): 105 return x 106 } 107 d := math.Cosh(2*real(x)) + math.Cos(2*imag(x)) 108 if d == 0 { 109 return Inf() 110 } 111 return complex(math.Sinh(2*real(x))/d, math.Sin(2*imag(x))/d) 112 } 113 114 // reducePi reduces the input argument x to the range (-Pi/2, Pi/2]. 115 // x must be greater than or equal to 0. For small arguments it 116 // uses Cody-Waite reduction in 3 float64 parts based on: 117 // "Elementary Function Evaluation: Algorithms and Implementation" 118 // Jean-Michel Muller, 1997. 119 // For very large arguments it uses Payne-Hanek range reduction based on: 120 // "ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit" 121 // K. C. Ng et al, March 24, 1992. 122 func reducePi(x float64) float64 { 123 // reduceThreshold is the maximum value of x where the reduction using 124 // Cody-Waite reduction still gives accurate results. This threshold 125 // is set by t*PIn being representable as a float64 without error 126 // where t is given by t = floor(x * (1 / Pi)) and PIn are the leading partial 127 // terms of Pi. Since the leading terms, PI1 and PI2 below, have 30 and 32 128 // trailing zero bits respectively, t should have less than 30 significant bits. 129 // t < 1<<30 -> floor(x*(1/Pi)+0.5) < 1<<30 -> x < (1<<30-1) * Pi - 0.5 130 // So, conservatively we can take x < 1<<30. 131 const reduceThreshold float64 = 1 << 30 132 if math.Abs(x) < reduceThreshold { 133 // Use Cody-Waite reduction in three parts. 134 const ( 135 // PI1, PI2 and PI3 comprise an extended precision value of PI 136 // such that PI ~= PI1 + PI2 + PI3. The parts are chosen so 137 // that PI1 and PI2 have an approximately equal number of trailing 138 // zero bits. This ensures that t*PI1 and t*PI2 are exact for 139 // large integer values of t. The full precision PI3 ensures the 140 // approximation of PI is accurate to 102 bits to handle cancellation 141 // during subtraction. 142 PI1 = 3.141592502593994 // 0x400921fb40000000 143 PI2 = 1.5099578831723193e-07 // 0x3e84442d00000000 144 PI3 = 1.0780605716316238e-14 // 0x3d08469898cc5170 145 ) 146 t := x / math.Pi 147 t += 0.5 148 t = float64(int64(t)) // int64(t) = the multiple 149 return ((x - t*PI1) - t*PI2) - t*PI3 150 } 151 // Must apply Payne-Hanek range reduction 152 const ( 153 mask = 0x7FF 154 shift = 64 - 11 - 1 155 bias = 1023 156 fracMask = 1<<shift - 1 157 ) 158 // Extract out the integer and exponent such that, 159 // x = ix * 2 ** exp. 160 ix := math.Float64bits(x) 161 exp := int(ix>>shift&mask) - bias - shift 162 ix &= fracMask 163 ix |= 1 << shift 164 165 // mPi is the binary digits of 1/Pi as a uint64 array, 166 // that is, 1/Pi = Sum mPi[i]*2^(-64*i). 167 // 19 64-bit digits give 1216 bits of precision 168 // to handle the largest possible float64 exponent. 169 var mPi = [...]uint64{ 170 0x0000000000000000, 171 0x517cc1b727220a94, 172 0xfe13abe8fa9a6ee0, 173 0x6db14acc9e21c820, 174 0xff28b1d5ef5de2b0, 175 0xdb92371d2126e970, 176 0x0324977504e8c90e, 177 0x7f0ef58e5894d39f, 178 0x74411afa975da242, 179 0x74ce38135a2fbf20, 180 0x9cc8eb1cc1a99cfa, 181 0x4e422fc5defc941d, 182 0x8ffc4bffef02cc07, 183 0xf79788c5ad05368f, 184 0xb69b3f6793e584db, 185 0xa7a31fb34f2ff516, 186 0xba93dd63f5f2f8bd, 187 0x9e839cfbc5294975, 188 0x35fdafd88fc6ae84, 189 0x2b0198237e3db5d5, 190 } 191 // Use the exponent to extract the 3 appropriate uint64 digits from mPi, 192 // B ~ (z0, z1, z2), such that the product leading digit has the exponent -64. 193 // Note, exp >= 50 since x >= reduceThreshold and exp < 971 for maximum float64. 194 digit, bitshift := uint(exp+64)/64, uint(exp+64)%64 195 z0 := (mPi[digit] << bitshift) | (mPi[digit+1] >> (64 - bitshift)) 196 z1 := (mPi[digit+1] << bitshift) | (mPi[digit+2] >> (64 - bitshift)) 197 z2 := (mPi[digit+2] << bitshift) | (mPi[digit+3] >> (64 - bitshift)) 198 // Multiply mantissa by the digits and extract the upper two digits (hi, lo). 199 z2hi, _ := bits.Mul64(z2, ix) 200 z1hi, z1lo := bits.Mul64(z1, ix) 201 z0lo := z0 * ix 202 lo, c := bits.Add64(z1lo, z2hi, 0) 203 hi, _ := bits.Add64(z0lo, z1hi, c) 204 // Find the magnitude of the fraction. 205 lz := uint(bits.LeadingZeros64(hi)) 206 e := uint64(bias - (lz + 1)) 207 // Clear implicit mantissa bit and shift into place. 208 hi = (hi << (lz + 1)) | (lo >> (64 - (lz + 1))) 209 hi >>= 64 - shift 210 // Include the exponent and convert to a float. 211 hi |= e << shift 212 x = math.Float64frombits(hi) 213 // map to (-Pi/2, Pi/2] 214 if x > 0.5 { 215 x-- 216 } 217 return math.Pi * x 218 } 219 220 // Taylor series expansion for cosh(2y) - cos(2x) 221 func tanSeries(z complex128) float64 { 222 const MACHEP = 1.0 / (1 << 53) 223 x := math.Abs(2 * real(z)) 224 y := math.Abs(2 * imag(z)) 225 x = reducePi(x) 226 x = x * x 227 y = y * y 228 x2 := 1.0 229 y2 := 1.0 230 f := 1.0 231 rn := 0.0 232 d := 0.0 233 for { 234 rn++ 235 f *= rn 236 rn++ 237 f *= rn 238 x2 *= x 239 y2 *= y 240 t := y2 + x2 241 t /= f 242 d += t 243 244 rn++ 245 f *= rn 246 rn++ 247 f *= rn 248 x2 *= x 249 y2 *= y 250 t = y2 - x2 251 t /= f 252 d += t 253 if !(math.Abs(t/d) > MACHEP) { 254 // Caution: Use ! and > instead of <= for correct behavior if t/d is NaN. 255 // See issue 17577. 256 break 257 } 258 } 259 return d 260 } 261 262 // Complex circular cotangent 263 // 264 // DESCRIPTION: 265 // 266 // If 267 // z = x + iy, 268 // 269 // then 270 // 271 // sin 2x - i sinh 2y 272 // w = --------------------. 273 // cosh 2y - cos 2x 274 // 275 // On the real axis, the denominator has zeros at even 276 // multiples of PI/2. Near these points it is evaluated 277 // by a Taylor series. 278 // 279 // ACCURACY: 280 // 281 // Relative error: 282 // arithmetic domain # trials peak rms 283 // DEC -10,+10 3000 6.5e-17 1.6e-17 284 // IEEE -10,+10 30000 9.2e-16 1.2e-16 285 // Also tested by ctan * ccot = 1 + i0. 286 287 // Cot returns the cotangent of x. 288 func Cot(x complex128) complex128 { 289 d := math.Cosh(2*imag(x)) - math.Cos(2*real(x)) 290 if math.Abs(d) < 0.25 { 291 d = tanSeries(x) 292 } 293 if d == 0 { 294 return Inf() 295 } 296 return complex(math.Sin(2*real(x))/d, -math.Sinh(2*imag(x))/d) 297 } 298