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 // moshier@na-net.ornl.gov
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