Source file src/pkg/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 "math"
     8	
     9	// The original C code, the long comment, and the constants
    10	// below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
    11	// The go code is a simplified version of the original C.
    12	//
    13	// Cephes Math Library Release 2.8:  June, 2000
    14	// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
    15	//
    16	// The readme file at http://netlib.sandia.gov/cephes/ says:
    17	//    Some software in this archive may be from the book _Methods and
    18	// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
    19	// International, 1989) or from the Cephes Mathematical Library, a
    20	// commercial product. In either event, it is copyrighted by the author.
    21	// What you see here may be used freely but it comes with no support or
    22	// guarantee.
    23	//
    24	//   The two known misprints in the book are repaired here in the
    25	// source listings for the gamma function and the incomplete beta
    26	// integral.
    27	//
    28	//   Stephen L. Moshier
    29	//   [email protected]
    30	
    31	// Complex circular tangent
    32	//
    33	// DESCRIPTION:
    34	//
    35	// If
    36	//     z = x + iy,
    37	//
    38	// then
    39	//
    40	//           sin 2x  +  i sinh 2y
    41	//     w  =  --------------------.
    42	//            cos 2x  +  cosh 2y
    43	//
    44	// On the real axis the denominator is zero at odd multiples
    45	// of PI/2.  The denominator is evaluated by its Taylor
    46	// series near these points.
    47	//
    48	// ctan(z) = -i ctanh(iz).
    49	//
    50	// ACCURACY:
    51	//
    52	//                      Relative error:
    53	// arithmetic   domain     # trials      peak         rms
    54	//    DEC       -10,+10      5200       7.1e-17     1.6e-17
    55	//    IEEE      -10,+10     30000       7.2e-16     1.2e-16
    56	// Also tested by ctan * ccot = 1 and catan(ctan(z))  =  z.
    57	
    58	// Tan returns the tangent of x.
    59	func Tan(x complex128) complex128 {
    60		d := math.Cos(2*real(x)) + math.Cosh(2*imag(x))
    61		if math.Abs(d) < 0.25 {
    62			d = tanSeries(x)
    63		}
    64		if d == 0 {
    65			return Inf()
    66		}
    67		return complex(math.Sin(2*real(x))/d, math.Sinh(2*imag(x))/d)
    68	}
    69	
    70	// Complex hyperbolic tangent
    71	//
    72	// DESCRIPTION:
    73	//
    74	// tanh z = (sinh 2x  +  i sin 2y) / (cosh 2x + cos 2y) .
    75	//
    76	// ACCURACY:
    77	//
    78	//                      Relative error:
    79	// arithmetic   domain     # trials      peak         rms
    80	//    IEEE      -10,+10     30000       1.7e-14     2.4e-16
    81	
    82	// Tanh returns the hyperbolic tangent of x.
    83	func Tanh(x complex128) complex128 {
    84		d := math.Cosh(2*real(x)) + math.Cos(2*imag(x))
    85		if d == 0 {
    86			return Inf()
    87		}
    88		return complex(math.Sinh(2*real(x))/d, math.Sin(2*imag(x))/d)
    89	}
    90	
    91	// Program to subtract nearest integer multiple of PI
    92	func reducePi(x float64) float64 {
    93		const (
    94			// extended precision value of PI:
    95			DP1 = 3.14159265160560607910E0   // ?? 0x400921fb54000000
    96			DP2 = 1.98418714791870343106E-9  // ?? 0x3e210b4610000000
    97			DP3 = 1.14423774522196636802E-17 // ?? 0x3c6a62633145c06e
    98		)
    99		t := x / math.Pi
   100		if t >= 0 {
   101			t += 0.5
   102		} else {
   103			t -= 0.5
   104		}
   105		t = float64(int64(t)) // int64(t) = the multiple
   106		return ((x - t*DP1) - t*DP2) - t*DP3
   107	}
   108	
   109	// Taylor series expansion for cosh(2y) - cos(2x)
   110	func tanSeries(z complex128) float64 {
   111		const MACHEP = 1.0 / (1 << 53)
   112		x := math.Abs(2 * real(z))
   113		y := math.Abs(2 * imag(z))
   114		x = reducePi(x)
   115		x = x * x
   116		y = y * y
   117		x2 := 1.0
   118		y2 := 1.0
   119		f := 1.0
   120		rn := 0.0
   121		d := 0.0
   122		for {
   123			rn += 1
   124			f *= rn
   125			rn += 1
   126			f *= rn
   127			x2 *= x
   128			y2 *= y
   129			t := y2 + x2
   130			t /= f
   131			d += t
   132	
   133			rn += 1
   134			f *= rn
   135			rn += 1
   136			f *= rn
   137			x2 *= x
   138			y2 *= y
   139			t = y2 - x2
   140			t /= f
   141			d += t
   142			if math.Abs(t/d) <= MACHEP {
   143				break
   144			}
   145		}
   146		return d
   147	}
   148	
   149	// Complex circular cotangent
   150	//
   151	// DESCRIPTION:
   152	//
   153	// If
   154	//     z = x + iy,
   155	//
   156	// then
   157	//
   158	//           sin 2x  -  i sinh 2y
   159	//     w  =  --------------------.
   160	//            cosh 2y  -  cos 2x
   161	//
   162	// On the real axis, the denominator has zeros at even
   163	// multiples of PI/2.  Near these points it is evaluated
   164	// by a Taylor series.
   165	//
   166	// ACCURACY:
   167	//
   168	//                      Relative error:
   169	// arithmetic   domain     # trials      peak         rms
   170	//    DEC       -10,+10      3000       6.5e-17     1.6e-17
   171	//    IEEE      -10,+10     30000       9.2e-16     1.2e-16
   172	// Also tested by ctan * ccot = 1 + i0.
   173	
   174	// Cot returns the cotangent of x.
   175	func Cot(x complex128) complex128 {
   176		d := math.Cosh(2*imag(x)) - math.Cos(2*real(x))
   177		if math.Abs(d) < 0.25 {
   178			d = tanSeries(x)
   179		}
   180		if d == 0 {
   181			return Inf()
   182		}
   183		return complex(math.Sin(2*real(x))/d, -math.Sinh(2*imag(x))/d)
   184	}