Source file src/pkg/math/tan.go

     1	// Copyright 2011 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 math
     6	
     7	/*
     8		Floating-point tangent.
     9	*/
    10	
    11	// The original C code, the long comment, and the constants
    12	// below were from http://netlib.sandia.gov/cephes/cmath/sin.c,
    13	// available from http://www.netlib.org/cephes/cmath.tgz.
    14	// The go code is a simplified version of the original C.
    15	//
    16	//      tan.c
    17	//
    18	//      Circular tangent
    19	//
    20	// SYNOPSIS:
    21	//
    22	// double x, y, tan();
    23	// y = tan( x );
    24	//
    25	// DESCRIPTION:
    26	//
    27	// Returns the circular tangent of the radian argument x.
    28	//
    29	// Range reduction is modulo pi/4.  A rational function
    30	//       x + x**3 P(x**2)/Q(x**2)
    31	// is employed in the basic interval [0, pi/4].
    32	//
    33	// ACCURACY:
    34	//                      Relative error:
    35	// arithmetic   domain     # trials      peak         rms
    36	//    DEC      +-1.07e9      44000      4.1e-17     1.0e-17
    37	//    IEEE     +-1.07e9      30000      2.9e-16     8.1e-17
    38	//
    39	// Partial loss of accuracy begins to occur at x = 2**30 = 1.074e9.  The loss
    40	// is not gradual, but jumps suddenly to about 1 part in 10e7.  Results may
    41	// be meaningless for x > 2**49 = 5.6e14.
    42	// [Accuracy loss statement from sin.go comments.]
    43	//
    44	// Cephes Math Library Release 2.8:  June, 2000
    45	// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
    46	//
    47	// The readme file at http://netlib.sandia.gov/cephes/ says:
    48	//    Some software in this archive may be from the book _Methods and
    49	// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
    50	// International, 1989) or from the Cephes Mathematical Library, a
    51	// commercial product. In either event, it is copyrighted by the author.
    52	// What you see here may be used freely but it comes with no support or
    53	// guarantee.
    54	//
    55	//   The two known misprints in the book are repaired here in the
    56	// source listings for the gamma function and the incomplete beta
    57	// integral.
    58	//
    59	//   Stephen L. Moshier
    60	//   [email protected]
    61	
    62	// tan coefficients
    63	var _tanP = [...]float64{
    64		-1.30936939181383777646E4, // 0xc0c992d8d24f3f38
    65		1.15351664838587416140E6,  // 0x413199eca5fc9ddd
    66		-1.79565251976484877988E7, // 0xc1711fead3299176
    67	}
    68	var _tanQ = [...]float64{
    69		1.00000000000000000000E0,
    70		1.36812963470692954678E4,  //0x40cab8a5eeb36572
    71		-1.32089234440210967447E6, //0xc13427bc582abc96
    72		2.50083801823357915839E7,  //0x4177d98fc2ead8ef
    73		-5.38695755929454629881E7, //0xc189afe03cbe5a31
    74	}
    75	
    76	// Tan returns the tangent of x.
    77	//
    78	// Special cases are:
    79	//	Tan(±0) = ±0
    80	//	Tan(±Inf) = NaN
    81	//	Tan(NaN) = NaN
    82	func Tan(x float64) float64
    83	
    84	func tan(x float64) float64 {
    85		const (
    86			PI4A = 7.85398125648498535156E-1                             // 0x3fe921fb40000000, Pi/4 split into three parts
    87			PI4B = 3.77489470793079817668E-8                             // 0x3e64442d00000000,
    88			PI4C = 2.69515142907905952645E-15                            // 0x3ce8469898cc5170,
    89			M4PI = 1.273239544735162542821171882678754627704620361328125 // 4/pi
    90		)
    91		// special cases
    92		switch {
    93		case x == 0 || IsNaN(x):
    94			return x // return ±0 || NaN()
    95		case IsInf(x, 0):
    96			return NaN()
    97		}
    98	
    99		// make argument positive but save the sign
   100		sign := false
   101		if x < 0 {
   102			x = -x
   103			sign = true
   104		}
   105	
   106		j := int64(x * M4PI) // integer part of x/(Pi/4), as integer for tests on the phase angle
   107		y := float64(j)      // integer part of x/(Pi/4), as float
   108	
   109		/* map zeros and singularities to origin */
   110		if j&1 == 1 {
   111			j += 1
   112			y += 1
   113		}
   114	
   115		z := ((x - y*PI4A) - y*PI4B) - y*PI4C
   116		zz := z * z
   117	
   118		if zz > 1e-14 {
   119			y = z + z*(zz*(((_tanP[0]*zz)+_tanP[1])*zz+_tanP[2])/((((zz+_tanQ[1])*zz+_tanQ[2])*zz+_tanQ[3])*zz+_tanQ[4]))
   120		} else {
   121			y = z
   122		}
   123		if j&2 == 2 {
   124			y = -1 / y
   125		}
   126		if sign {
   127			y = -y
   128		}
   129		return y
   130	}