libgo/go/math/cmplx/asin.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 //   moshier@na-net.ornl.gov
  30
  31 // Complex circular arc sine
  32 //
  33 // DESCRIPTION:
  34 //
  35 // Inverse complex sine:
  36 //                               2
  37 // w = -i clog( iz + csqrt( 1 - z ) ).
  38 //
  39 // casin(z) = -i casinh(iz)
  40 //
  41 // ACCURACY:
  42 //
  43 //                      Relative error:
  44 // arithmetic   domain     # trials      peak         rms
  45 //    DEC       -10,+10     10100       2.1e-15     3.4e-16
  46 //    IEEE      -10,+10     30000       2.2e-14     2.7e-15
  47 // Larger relative error can be observed for z near zero.
  48 // Also tested by csin(casin(z)) = z.
  49
  50 // Asin returns the inverse sine of x.
  51 func Asin(x complex128) complex128 {
  52         switch re, im := real(x), imag(x); {
  53         case im == 0 && math.Abs(re) <= 1:
  54                 return complex(math.Asin(re), im)
  55         case re == 0 && math.Abs(im) <= 1:
  56                 return complex(re, math.Asinh(im))
  57         case math.IsNaN(im):
  58                 switch {
  59                 case re == 0:
  60                         return complex(re, math.NaN())
  61                 case math.IsInf(re, 0):
  62                         return complex(math.NaN(), re)
  63                 default:
  64                         return NaN()
  65                 }
  66         case math.IsInf(im, 0):
  67                 switch {
  68                 case math.IsNaN(re):
  69                         return x
  70                 case math.IsInf(re, 0):
  71                         return complex(math.Copysign(math.Pi/4, re), im)
  72                 default:
  73                         return complex(math.Copysign(0, re), im)
  74                 }
  75         case math.IsInf(re, 0):
  76                 return complex(math.Copysign(math.Pi/2, re), math.Copysign(re, im))
  77         }
  78         ct := complex(-imag(x), real(x)) // i * x
  79         xx := x * x
  80         x1 := complex(1-real(xx), -imag(xx)) // 1 - x*x
  81         x2 := Sqrt(x1)                       // x2 = sqrt(1 - x*x)
  82         w := Log(ct + x2)
  83         return complex(imag(w), -real(w)) // -i * w
  84 }
  85
  86 // Asinh returns the inverse hyperbolic sine of x.
  87 func Asinh(x complex128) complex128 {
  88         switch re, im := real(x), imag(x); {
  89         case im == 0 && math.Abs(re) <= 1:
  90                 return complex(math.Asinh(re), im)
  91         case re == 0 && math.Abs(im) <= 1:
  92                 return complex(re, math.Asin(im))
  93         case math.IsInf(re, 0):
  94                 switch {
  95                 case math.IsInf(im, 0):
  96                         return complex(re, math.Copysign(math.Pi/4, im))
  97                 case math.IsNaN(im):
  98                         return x
  99                 default:
 100                         return complex(re, math.Copysign(0.0, im))
 101                 }
 102         case math.IsNaN(re):
 103                 switch {
 104                 case im == 0:
 105                         return x
 106                 case math.IsInf(im, 0):
 107                         return complex(im, re)
 108                 default:
 109                         return NaN()
 110                 }
 111         case math.IsInf(im, 0):
 112                 return complex(math.Copysign(im, re), math.Copysign(math.Pi/2, im))
 113         }
 114         xx := x * x
 115         x1 := complex(1+real(xx), imag(xx)) // 1 + x*x
 116         return Log(x + Sqrt(x1))            // log(x + sqrt(1 + x*x))
 117 }
 118
 119 // Complex circular arc cosine
 120 //
 121 // DESCRIPTION:
 122 //
 123 // w = arccos z  =  PI/2 - arcsin z.
 124 //
 125 // ACCURACY:
 126 //
 127 //                      Relative error:
 128 // arithmetic   domain     # trials      peak         rms
 129 //    DEC       -10,+10      5200      1.6e-15      2.8e-16
 130 //    IEEE      -10,+10     30000      1.8e-14      2.2e-15
 131
 132 // Acos returns the inverse cosine of x.
 133 func Acos(x complex128) complex128 {
 134         w := Asin(x)
 135         return complex(math.Pi/2-real(w), -imag(w))
 136 }
 137
 138 // Acosh returns the inverse hyperbolic cosine of x.
 139 func Acosh(x complex128) complex128 {
 140         if x == 0 {
 141                 return complex(0, math.Copysign(math.Pi/2, imag(x)))
 142         }
 143         w := Acos(x)
 144         if imag(w) <= 0 {
 145                 return complex(-imag(w), real(w)) // i * w
 146         }
 147         return complex(imag(w), -real(w)) // -i * w
 148 }
 149
 150 // Complex circular arc tangent
 151 //
 152 // DESCRIPTION:
 153 //
 154 // If
 155 //     z = x + iy,
 156 //
 157 // then
 158 //          1       (    2x     )
 159 // Re w  =  - arctan(-----------)  +  k PI
 160 //          2       (     2    2)
 161 //                  (1 - x  - y )
 162 //
 163 //               ( 2         2)
 164 //          1    (x  +  (y+1) )
 165 // Im w  =  - log(------------)
 166 //          4    ( 2         2)
 167 //               (x  +  (y-1) )
 168 //
 169 // Where k is an arbitrary integer.
 170 //
 171 // catan(z) = -i catanh(iz).
 172 //
 173 // ACCURACY:
 174 //
 175 //                      Relative error:
 176 // arithmetic   domain     # trials      peak         rms
 177 //    DEC       -10,+10      5900       1.3e-16     7.8e-18
 178 //    IEEE      -10,+10     30000       2.3e-15     8.5e-17
 179 // The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,
 180 // had peak relative error 1.5e-16, rms relative error
 181 // 2.9e-17.  See also clog().
 182
 183 // Atan returns the inverse tangent of x.
 184 func Atan(x complex128) complex128 {
 185         switch re, im := real(x), imag(x); {
 186         case im == 0:
 187                 return complex(math.Atan(re), im)
 188         case re == 0 && math.Abs(im) <= 1:
 189                 return complex(re, math.Atanh(im))
 190         case math.IsInf(im, 0) || math.IsInf(re, 0):
 191                 if math.IsNaN(re) {
 192                         return complex(math.NaN(), math.Copysign(0, im))
 193                 }
 194                 return complex(math.Copysign(math.Pi/2, re), math.Copysign(0, im))
 195         case math.IsNaN(re) || math.IsNaN(im):
 196                 return NaN()
 197         }
 198         x2 := real(x) * real(x)
 199         a := 1 - x2 - imag(x)*imag(x)
 200         if a == 0 {
 201                 return NaN()
 202         }
 203         t := 0.5 * math.Atan2(2*real(x), a)
 204         w := reducePi(t)
 205
 206         t = imag(x) - 1
 207         b := x2 + t*t
 208         if b == 0 {
 209                 return NaN()
 210         }
 211         t = imag(x) + 1
 212         c := (x2 + t*t) / b
 213         return complex(w, 0.25*math.Log(c))
 214 }
 215
 216 // Atanh returns the inverse hyperbolic tangent of x.
 217 func Atanh(x complex128) complex128 {
 218         z := complex(-imag(x), real(x)) // z = i * x
 219         z = Atan(z)
 220         return complex(imag(z), -real(z)) // z = -i * z
 221 }