libgo/go/cmath/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 cmath
   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 //   moshier@na-net.ornl.gov
  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.Fabs(d) < 0.25 {
  62                 d = tanSeries(x)
  63         }
  64         if d == 0 {
  65                 return Inf()
  66         }
  67         return cmplx(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 cmplx(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.Fabs(2 * real(z))
 113         y := math.Fabs(2 * imag(z))
 114         x = reducePi(x)
 115         x = x * x
 116         y = y * y
 117         x2 := float64(1)
 118         y2 := float64(1)
 119         f := float64(1)
 120         rn := float64(0)
 121         d := float64(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.Fabs(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.Fabs(d) < 0.25 {
 178                 d = tanSeries(x)
 179         }
 180         if d == 0 {
 181                 return Inf()
 182         }
 183         return cmplx(math.Sin(2*real(x))/d, -math.Sinh(2*imag(x))/d)
 184 }