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
     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  

View as plain text