...

Source file src/math/cmplx/tan.go

Documentation: math/cmplx

		 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  

View as plain text