...

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	//   moshier@na-net.ornl.gov
    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 the radian argument 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		)
    90		// special cases
    91		switch {
    92		case x == 0 || IsNaN(x):
    93			return x // return ±0 || NaN()
    94		case IsInf(x, 0):
    95			return NaN()
    96		}
    97	
    98		// make argument positive but save the sign
    99		sign := false
   100		if x < 0 {
   101			x = -x
   102			sign = true
   103		}
   104		var j uint64
   105		var y, z float64
   106		if x >= reduceThreshold {
   107			j, z = trigReduce(x)
   108		} else {
   109			j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
   110			y = float64(j)           // integer part of x/(Pi/4), as float
   111	
   112			/* map zeros and singularities to origin */
   113			if j&1 == 1 {
   114				j++
   115				y++
   116			}
   117	
   118			z = ((x - y*PI4A) - y*PI4B) - y*PI4C
   119		}
   120		zz := z * z
   121	
   122		if zz > 1e-14 {
   123			y = z + z*(zz*(((_tanP[0]*zz)+_tanP[1])*zz+_tanP[2])/((((zz+_tanQ[1])*zz+_tanQ[2])*zz+_tanQ[3])*zz+_tanQ[4]))
   124		} else {
   125			y = z
   126		}
   127		if j&2 == 2 {
   128			y = -1 / y
   129		}
   130		if sign {
   131			y = -y
   132		}
   133		return y
   134	}
   135	

View as plain text