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 5package cmplx 6 7import "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. 53func 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. 86func 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. 124func 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. 156func 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 176func 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