tan.go 5.03 KB
Newer Older
Charles L. Dorian's avatar
Charles L. Dorian committed
1 2 3 4
// Copyright 2010 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

5
package cmplx
Charles L. Dorian's avatar
Charles L. Dorian committed
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

import "math"

// The original C code, the long comment, and the constants
// below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
// The go code is a simplified version of the original C.
//
// Cephes Math Library Release 2.8:  June, 2000
// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
//
// The readme file at http://netlib.sandia.gov/cephes/ says:
//    Some software in this archive may be from the book _Methods and
// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
// International, 1989) or from the Cephes Mathematical Library, a
// commercial product. In either event, it is copyrighted by the author.
// What you see here may be used freely but it comes with no support or
// guarantee.
//
//   The two known misprints in the book are repaired here in the
// source listings for the gamma function and the incomplete beta
// integral.
//
//   Stephen L. Moshier
//   moshier@na-net.ornl.gov

// Complex circular tangent
//
// DESCRIPTION:
//
// If
//     z = x + iy,
//
// then
//
//           sin 2x  +  i sinh 2y
//     w  =  --------------------.
//            cos 2x  +  cosh 2y
//
// On the real axis the denominator is zero at odd multiples
// of PI/2.  The denominator is evaluated by its Taylor
// series near these points.
//
// ctan(z) = -i ctanh(iz).
//
// ACCURACY:
//
//                      Relative error:
// arithmetic   domain     # trials      peak         rms
//    DEC       -10,+10      5200       7.1e-17     1.6e-17
//    IEEE      -10,+10     30000       7.2e-16     1.2e-16
// Also tested by ctan * ccot = 1 and catan(ctan(z))  =  z.

// Tan returns the tangent of x.
func Tan(x complex128) complex128 {
60 61 62 63 64 65 66 67 68 69
	switch re, im := real(x), imag(x); {
	case math.IsInf(im, 0):
		switch {
		case math.IsInf(re, 0) || math.IsNaN(re):
			return complex(math.Copysign(0, re), math.Copysign(1, im))
		}
		return complex(math.Copysign(0, math.Sin(2*re)), math.Copysign(1, im))
	case re == 0 && math.IsNaN(im):
		return x
	}
Charles L. Dorian's avatar
Charles L. Dorian committed
70
	d := math.Cos(2*real(x)) + math.Cosh(2*imag(x))
71
	if math.Abs(d) < 0.25 {
Charles L. Dorian's avatar
Charles L. Dorian committed
72 73 74 75 76
		d = tanSeries(x)
	}
	if d == 0 {
		return Inf()
	}
77
	return complex(math.Sin(2*real(x))/d, math.Sinh(2*imag(x))/d)
Charles L. Dorian's avatar
Charles L. Dorian committed
78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
}

// Complex hyperbolic tangent
//
// DESCRIPTION:
//
// tanh z = (sinh 2x  +  i sin 2y) / (cosh 2x + cos 2y) .
//
// ACCURACY:
//
//                      Relative error:
// arithmetic   domain     # trials      peak         rms
//    IEEE      -10,+10     30000       1.7e-14     2.4e-16

// Tanh returns the hyperbolic tangent of x.
func Tanh(x complex128) complex128 {
94 95 96 97 98 99 100 101 102 103
	switch re, im := real(x), imag(x); {
	case math.IsInf(re, 0):
		switch {
		case math.IsInf(im, 0) || math.IsNaN(im):
			return complex(math.Copysign(1, re), math.Copysign(0, im))
		}
		return complex(math.Copysign(1, re), math.Copysign(0, math.Sin(2*im)))
	case im == 0 && math.IsNaN(re):
		return x
	}
Charles L. Dorian's avatar
Charles L. Dorian committed
104 105 106 107
	d := math.Cosh(2*real(x)) + math.Cos(2*imag(x))
	if d == 0 {
		return Inf()
	}
108
	return complex(math.Sinh(2*real(x))/d, math.Sin(2*imag(x))/d)
Charles L. Dorian's avatar
Charles L. Dorian committed
109 110 111 112 113 114
}

// Program to subtract nearest integer multiple of PI
func reducePi(x float64) float64 {
	const (
		// extended precision value of PI:
Robert Griesemer's avatar
Robert Griesemer committed
115 116 117
		DP1 = 3.14159265160560607910e0   // ?? 0x400921fb54000000
		DP2 = 1.98418714791870343106e-9  // ?? 0x3e210b4610000000
		DP3 = 1.14423774522196636802e-17 // ?? 0x3c6a62633145c06e
Charles L. Dorian's avatar
Charles L. Dorian committed
118 119 120 121 122 123 124 125 126 127 128 129 130 131
	)
	t := x / math.Pi
	if t >= 0 {
		t += 0.5
	} else {
		t -= 0.5
	}
	t = float64(int64(t)) // int64(t) = the multiple
	return ((x - t*DP1) - t*DP2) - t*DP3
}

// Taylor series expansion for cosh(2y) - cos(2x)
func tanSeries(z complex128) float64 {
	const MACHEP = 1.0 / (1 << 53)
132 133
	x := math.Abs(2 * real(z))
	y := math.Abs(2 * imag(z))
Charles L. Dorian's avatar
Charles L. Dorian committed
134 135 136
	x = reducePi(x)
	x = x * x
	y = y * y
137 138 139 140 141
	x2 := 1.0
	y2 := 1.0
	f := 1.0
	rn := 0.0
	d := 0.0
Charles L. Dorian's avatar
Charles L. Dorian committed
142
	for {
143
		rn++
Charles L. Dorian's avatar
Charles L. Dorian committed
144
		f *= rn
145
		rn++
Charles L. Dorian's avatar
Charles L. Dorian committed
146 147 148 149 150 151 152
		f *= rn
		x2 *= x
		y2 *= y
		t := y2 + x2
		t /= f
		d += t

153
		rn++
Charles L. Dorian's avatar
Charles L. Dorian committed
154
		f *= rn
155
		rn++
Charles L. Dorian's avatar
Charles L. Dorian committed
156 157 158 159 160 161
		f *= rn
		x2 *= x
		y2 *= y
		t = y2 - x2
		t /= f
		d += t
162 163 164
		if !(math.Abs(t/d) > MACHEP) {
			// Caution: Use ! and > instead of <= for correct behavior if t/d is NaN.
			// See issue 17577.
Charles L. Dorian's avatar
Charles L. Dorian committed
165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
			break
		}
	}
	return d
}

// Complex circular cotangent
//
// DESCRIPTION:
//
// If
//     z = x + iy,
//
// then
//
//           sin 2x  -  i sinh 2y
//     w  =  --------------------.
//            cosh 2y  -  cos 2x
//
// On the real axis, the denominator has zeros at even
// multiples of PI/2.  Near these points it is evaluated
// by a Taylor series.
//
// ACCURACY:
//
//                      Relative error:
// arithmetic   domain     # trials      peak         rms
//    DEC       -10,+10      3000       6.5e-17     1.6e-17
//    IEEE      -10,+10     30000       9.2e-16     1.2e-16
// Also tested by ctan * ccot = 1 + i0.

// Cot returns the cotangent of x.
func Cot(x complex128) complex128 {
	d := math.Cosh(2*imag(x)) - math.Cos(2*real(x))
199
	if math.Abs(d) < 0.25 {
Charles L. Dorian's avatar
Charles L. Dorian committed
200 201 202 203 204
		d = tanSeries(x)
	}
	if d == 0 {
		return Inf()
	}
205
	return complex(math.Sin(2*real(x))/d, -math.Sinh(2*imag(x))/d)
Charles L. Dorian's avatar
Charles L. Dorian committed
206
}