# Source file src/math/cmplx/tan.go

```     1  // Copyright 2010 The Go Authors. All rights reserved.
2  // Use of this source code is governed by a BSD-style
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 (
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
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,
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.
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
```

View as plain text