sin.mx raw

   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 sine
  32  //
  33  // DESCRIPTION:
  34  //
  35  // If
  36  //     z = x + iy,
  37  //
  38  // then
  39  //
  40  //     w = sin x  cosh y  +  i cos x sinh y.
  41  //
  42  // csin(z) = -i csinh(iz).
  43  //
  44  // ACCURACY:
  45  //
  46  //                      Relative error:
  47  // arithmetic   domain     # trials      peak         rms
  48  //    DEC       -10,+10      8400       5.3e-17     1.3e-17
  49  //    IEEE      -10,+10     30000       3.8e-16     1.0e-16
  50  // Also tested by csin(casin(z)) = z.
  51  
  52  // Sin returns the sine of x.
  53  func Sin(x complex128) complex128 {
  54  	switch re, im := real(x), imag(x); {
  55  	case im == 0 && (math.IsInf(re, 0) || math.IsNaN(re)):
  56  		return complex(math.NaN(), im)
  57  	case math.IsInf(im, 0):
  58  		switch {
  59  		case re == 0:
  60  			return x
  61  		case math.IsInf(re, 0) || math.IsNaN(re):
  62  			return complex(math.NaN(), im)
  63  		}
  64  	case re == 0 && math.IsNaN(im):
  65  		return x
  66  	}
  67  	s, c := math.Sincos(real(x))
  68  	sh, ch := sinhcosh(imag(x))
  69  	return complex(s*ch, c*sh)
  70  }
  71  
  72  // Complex hyperbolic sine
  73  //
  74  // DESCRIPTION:
  75  //
  76  // csinh z = (cexp(z) - cexp(-z))/2
  77  //         = sinh x * cos y  +  i cosh x * sin y .
  78  //
  79  // ACCURACY:
  80  //
  81  //                      Relative error:
  82  // arithmetic   domain     # trials      peak         rms
  83  //    IEEE      -10,+10     30000       3.1e-16     8.2e-17
  84  
  85  // Sinh returns the hyperbolic sine of x.
  86  func Sinh(x complex128) complex128 {
  87  	switch re, im := real(x), imag(x); {
  88  	case re == 0 && (math.IsInf(im, 0) || math.IsNaN(im)):
  89  		return complex(re, math.NaN())
  90  	case math.IsInf(re, 0):
  91  		switch {
  92  		case im == 0:
  93  			return complex(re, im)
  94  		case math.IsInf(im, 0) || math.IsNaN(im):
  95  			return complex(re, math.NaN())
  96  		}
  97  	case im == 0 && math.IsNaN(re):
  98  		return complex(math.NaN(), im)
  99  	}
 100  	s, c := math.Sincos(imag(x))
 101  	sh, ch := sinhcosh(real(x))
 102  	return complex(c*sh, s*ch)
 103  }
 104  
 105  // Complex circular cosine
 106  //
 107  // DESCRIPTION:
 108  //
 109  // If
 110  //     z = x + iy,
 111  //
 112  // then
 113  //
 114  //     w = cos x  cosh y  -  i sin x sinh y.
 115  //
 116  // ACCURACY:
 117  //
 118  //                      Relative error:
 119  // arithmetic   domain     # trials      peak         rms
 120  //    DEC       -10,+10      8400       4.5e-17     1.3e-17
 121  //    IEEE      -10,+10     30000       3.8e-16     1.0e-16
 122  
 123  // Cos returns the cosine of x.
 124  func Cos(x complex128) complex128 {
 125  	switch re, im := real(x), imag(x); {
 126  	case im == 0 && (math.IsInf(re, 0) || math.IsNaN(re)):
 127  		return complex(math.NaN(), -im*math.Copysign(0, re))
 128  	case math.IsInf(im, 0):
 129  		switch {
 130  		case re == 0:
 131  			return complex(math.Inf(1), -re*math.Copysign(0, im))
 132  		case math.IsInf(re, 0) || math.IsNaN(re):
 133  			return complex(math.Inf(1), math.NaN())
 134  		}
 135  	case re == 0 && math.IsNaN(im):
 136  		return complex(math.NaN(), 0)
 137  	}
 138  	s, c := math.Sincos(real(x))
 139  	sh, ch := sinhcosh(imag(x))
 140  	return complex(c*ch, -s*sh)
 141  }
 142  
 143  // Complex hyperbolic cosine
 144  //
 145  // DESCRIPTION:
 146  //
 147  // ccosh(z) = cosh x  cos y + i sinh x sin y .
 148  //
 149  // ACCURACY:
 150  //
 151  //                      Relative error:
 152  // arithmetic   domain     # trials      peak         rms
 153  //    IEEE      -10,+10     30000       2.9e-16     8.1e-17
 154  
 155  // Cosh returns the hyperbolic cosine of x.
 156  func Cosh(x complex128) complex128 {
 157  	switch re, im := real(x), imag(x); {
 158  	case re == 0 && (math.IsInf(im, 0) || math.IsNaN(im)):
 159  		return complex(math.NaN(), re*math.Copysign(0, im))
 160  	case math.IsInf(re, 0):
 161  		switch {
 162  		case im == 0:
 163  			return complex(math.Inf(1), im*math.Copysign(0, re))
 164  		case math.IsInf(im, 0) || math.IsNaN(im):
 165  			return complex(math.Inf(1), math.NaN())
 166  		}
 167  	case im == 0 && math.IsNaN(re):
 168  		return complex(math.NaN(), im)
 169  	}
 170  	s, c := math.Sincos(imag(x))
 171  	sh, ch := sinhcosh(real(x))
 172  	return complex(c*ch, s*sh)
 173  }
 174  
 175  // calculate sinh and cosh.
 176  func sinhcosh(x float64) (sh, ch float64) {
 177  	if math.Abs(x) <= 0.5 {
 178  		return math.Sinh(x), math.Cosh(x)
 179  	}
 180  	e := math.Exp(x)
 181  	ei := 0.5 / e
 182  	e *= 0.5
 183  	return e - ei, e + ei
 184  }
 185