cmath: new package

Complex math function package. Still needs more special case checking.

R=rsc
CC=golang-dev
https://golang.org/cl/874041

src/pkg/Makefile

@@ -26,6 +26,7 @@ DIRS=\
 	bignum\
 	bufio\
 	bytes\
+	cmath\
 	compress/flate\
 	compress/gzip\
 	compress/zlib\

src/pkg/cmath/Makefile

+# Copyright 2010 The Go Authors. All rights reserved.
+# Use of this source code is governed by a BSD-style
+# license that can be found in the LICENSE file.
+
+include ../../Make.$(GOARCH)
+
+TARG=cmath
+
+GOFILES=\
+	abs.go\
+	asin.go\
+	conj.go\
+	exp.go\
+	isinf.go\
+	isnan.go\
+	log.go\
+	phase.go\
+	polar.go\
+	pow.go\
+	rect.go\
+	sin.go\
+	sqrt.go\
+	tan.go\
+
+include ../../Make.pkg

src/pkg/cmath/abs.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// Abs returns the absolute value (also called the modulus) of x.
+func Abs(x complex128) float64 { return math.Hypot(real(x), imag(x)) }

src/pkg/cmath/asin.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// The original C code, the long comment, and the constants
+// below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
+// The go code is a simplified version of the original C.
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   moshier@na-net.ornl.gov
+
+// Complex circular arc sine
+//
+// DESCRIPTION:
+//
+// Inverse complex sine:
+//                               2
+// w = -i clog( iz + csqrt( 1 - z ) ).
+//
+// casin(z) = -i casinh(iz)
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    DEC       -10,+10     10100       2.1e-15     3.4e-16
+//    IEEE      -10,+10     30000       2.2e-14     2.7e-15
+// Larger relative error can be observed for z near zero.
+// Also tested by csin(casin(z)) = z.
+
+// Asin returns the inverse sine of x.
+func Asin(x complex128) complex128 {
+	if imag(x) == 0 {
+		if math.Fabs(real(x)) > 1 {
+			return cmplx(math.Pi/2, 0) // DOMAIN error
+		}
+		return cmplx(math.Asin(real(x)), 0)
+	}
+	ct := cmplx(-imag(x), real(x)) // i * x
+	xx := x * x
+	x1 := cmplx(1-real(xx), -imag(xx)) // 1 - x*x
+	x2 := Sqrt(x1)                     // x2 = sqrt(1 - x*x)
+	w := Log(ct + x2)
+	return cmplx(imag(w), -real(w)) // -i * w
+}
+
+// Asinh returns the inverse hyperbolic sine of x.
+func Asinh(x complex128) complex128 {
+	// TODO check range
+	if imag(x) == 0 {
+		if math.Fabs(real(x)) > 1 {
+			return cmplx(math.Pi/2, 0) // DOMAIN error
+		}
+		return cmplx(math.Asinh(real(x)), 0)
+	}
+	xx := x * x
+	x1 := cmplx(1+real(xx), imag(xx)) // 1 + x*x
+	return Log(x + Sqrt(x1))          // log(x + sqrt(1 + x*x))
+}
+
+// Complex circular arc cosine
+//
+// DESCRIPTION:
+//
+// w = arccos z  =  PI/2 - arcsin z.
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    DEC       -10,+10      5200      1.6e-15      2.8e-16
+//    IEEE      -10,+10     30000      1.8e-14      2.2e-15
+
+// Acos returns the inverse cosine of x.
+func Acos(x complex128) complex128 {
+	w := Asin(x)
+	return cmplx(math.Pi/2-real(w), -imag(w))
+}
+
+// Acosh returns the inverse hyperbolic cosine of x.
+func Acosh(x complex128) complex128 {
+	w := Acos(x)
+	if imag(w) <= 0 {
+		return cmplx(-imag(w), real(w)) // i * w
+	}
+	return cmplx(imag(w), -real(w)) // -i * w
+}
+
+// Complex circular arc tangent
+//
+// DESCRIPTION:
+//
+// If
+//     z = x + iy,
+//
+// then
+//          1       (    2x     )
+// Re w  =  - arctan(-----------)  +  k PI
+//          2       (     2    2)
+//                  (1 - x  - y )
+//
+//               ( 2         2)
+//          1    (x  +  (y+1) )
+// Im w  =  - log(------------)
+//          4    ( 2         2)
+//               (x  +  (y-1) )
+//
+// Where k is an arbitrary integer.
+//
+// catan(z) = -i catanh(iz).
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    DEC       -10,+10      5900       1.3e-16     7.8e-18
+//    IEEE      -10,+10     30000       2.3e-15     8.5e-17
+// The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,
+// had peak relative error 1.5e-16, rms relative error
+// 2.9e-17.  See also clog().
+
+// Atan returns the inverse tangent of x.
+func Atan(x complex128) complex128 {
+	if real(x) == 0 && imag(x) > 1 {
+		return NaN()
+	}
+
+	x2 := real(x) * real(x)
+	a := 1 - x2 - imag(x)*imag(x)
+	if a == 0 {
+		return NaN()
+	}
+	t := 0.5 * math.Atan2(2*real(x), a)
+	w := reducePi(t)
+
+	t = imag(x) - 1
+	b := x2 + t*t
+	if b == 0 {
+		return NaN()
+	}
+	t = imag(x) + 1
+	c := (x2 + t*t) / b
+	return cmplx(w, 0.25*math.Log(c))
+}
+
+// Atanh returns the inverse hyperbolic tangent of x.
+func Atanh(x complex128) complex128 {
+	z := cmplx(-imag(x), real(x)) // z = i * x
+	z = Atan(z)
+	return cmplx(imag(z), -real(z)) // z = -i * z
+}

src/pkg/cmath/conj.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+// Conj returns the complex conjugate of x.
+func Conj(x complex128) complex128 { return cmplx(real(x), -imag(x)) }

src/pkg/cmath/exp.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// The original C code, the long comment, and the constants
+// below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
+// The go code is a simplified version of the original C.
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   moshier@na-net.ornl.gov
+
+// Complex exponential function
+//
+// DESCRIPTION:
+//
+// Returns the complex exponential of the complex argument z.
+//
+// If
+//     z = x + iy,
+//     r = exp(x),
+// then
+//     w = r cos y + i r sin y.
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    DEC       -10,+10      8700       3.7e-17     1.1e-17
+//    IEEE      -10,+10     30000       3.0e-16     8.7e-17
+
+// Exp returns e^x, the base-e exponential of x.
+func Exp(x complex128) complex128 {
+	r := math.Exp(real(x))
+	s, c := math.Sincos(imag(x))
+	return cmplx(r*c, r*s)
+}

src/pkg/cmath/isinf.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// IsInf returns true if either real(x) or imag(x) is an infinity.
+func IsInf(x complex128) bool {
+	if math.IsInf(real(x), 0) || math.IsInf(imag(x), 0) {
+		return true
+	}
+	return false
+}
+
+// Inf returns a complex infinity, cmplx(+Inf, +Inf).
+func Inf() complex128 {
+	inf := math.Inf(1)
+	return cmplx(inf, inf)
+}

src/pkg/cmath/isnan.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// IsNaN returns true if either real(x) or imag(x) is NaN.
+func IsNaN(x complex128) bool {
+	if math.IsNaN(real(x)) || math.IsNaN(imag(x)) {
+		return true
+	}
+	return false
+}
+
+// NaN returns a complex ``not-a-number'' value.
+func NaN() complex128 {
+	nan := math.NaN()
+	return cmplx(nan, nan)
+}

src/pkg/cmath/log.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// The original C code, the long comment, and the constants
+// below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
+// The go code is a simplified version of the original C.
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   moshier@na-net.ornl.gov
+
+// Complex natural logarithm
+//
+// DESCRIPTION:
+//
+// Returns complex logarithm to the base e (2.718...) of
+// the complex argument z.
+//
+// If
+//       z = x + iy, r = sqrt( x**2 + y**2 ),
+// then
+//       w = log(r) + i arctan(y/x).
+//
+// The arctangent ranges from -PI to +PI.
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    DEC       -10,+10      7000       8.5e-17     1.9e-17
+//    IEEE      -10,+10     30000       5.0e-15     1.1e-16
+//
+// Larger relative error can be observed for z near 1 +i0.
+// In IEEE arithmetic the peak absolute error is 5.2e-16, rms
+// absolute error 1.0e-16.
+
+// Log returns the natural logarithm of x.
+func Log(x complex128) complex128 {
+	return cmplx(math.Log(Abs(x)), Phase(x))
+}
+
+// Log10 returns the decimal logarithm of x.
+func Log10(x complex128) complex128 {
+	return math.Log10E * Log(x)
+}

src/pkg/cmath/phase.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// Phase returns the phase (also called the argument) of x.
+// The returned value is in the range (-Pi, Pi].
+func Phase(x complex128) float64 { return math.Atan2(imag(x), real(x)) }

src/pkg/cmath/polar.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+// Polar returns the absolute value r and phase θ of x,
+// such that x = r * e^θi.
+// The phase is in the range (-Pi, Pi].
+func Polar(x complex128) (r, θ float64) {
+	return Abs(x), Phase(x)
+}

src/pkg/cmath/pow.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// The original C code, the long comment, and the constants
+// below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
+// The go code is a simplified version of the original C.
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   moshier@na-net.ornl.gov
+
+// Complex power function
+//
+// DESCRIPTION:
+//
+// Raises complex A to the complex Zth power.
+// Definition is per AMS55 # 4.2.8,
+// analytically equivalent to cpow(a,z) = cexp(z clog(a)).
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    IEEE      -10,+10     30000       9.4e-15     1.5e-15
+
+// Pow returns x^y, the base-x exponential of y.
+func Pow(x, y complex128) complex128 {
+	modulus := Abs(x)
+	if modulus == 0 {
+		return cmplx(0, 0)
+	}
+	r := math.Pow(modulus, real(y))
+	arg := Phase(x)
+	theta := real(y) * arg
+	if imag(y) != 0 {
+		r *= math.Exp(-imag(y) * arg)
+		theta += imag(y) * math.Log(modulus)
+	}
+	s, c := math.Sincos(theta)
+	return cmplx(r*c, r*s)
+}

src/pkg/cmath/rect.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// Rect returns the complex number x with polar coordinates r, θ.
+func Rect(r, θ float64) complex128 {
+	s, c := math.Sincos(θ)
+	return cmplx(r*c, r*s)
+}

src/pkg/cmath/sin.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// The original C code, the long comment, and the constants
+// below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
+// The go code is a simplified version of the original C.
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   moshier@na-net.ornl.gov
+
+// Complex circular sine
+//
+// DESCRIPTION:
+//
+// If
+//     z = x + iy,
+//
+// then
+//
+//     w = sin x  cosh y  +  i cos x sinh y.
+//
+// csin(z) = -i csinh(iz).
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    DEC       -10,+10      8400       5.3e-17     1.3e-17
+//    IEEE      -10,+10     30000       3.8e-16     1.0e-16
+// Also tested by csin(casin(z)) = z.
+
+// Sin returns the sine of x.
+func Sin(x complex128) complex128 {
+	s, c := math.Sincos(real(x))
+	sh, ch := sinhcosh(imag(x))
+	return cmplx(s*ch, c*sh)
+}
+
+// Complex hyperbolic sine
+//
+// DESCRIPTION:
+//
+// csinh z = (cexp(z) - cexp(-z))/2
+//         = sinh x * cos y  +  i cosh x * sin y .
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    IEEE      -10,+10     30000       3.1e-16     8.2e-17
+
+// Sinh returns the hyperbolic sine of x.
+func Sinh(x complex128) complex128 {
+	s, c := math.Sincos(imag(x))
+	sh, ch := sinhcosh(real(x))
+	return cmplx(c*sh, s*ch)
+}
+
+// Complex circular cosine
+//
+// DESCRIPTION:
+//
+// If
+//     z = x + iy,
+//
+// then
+//
+//     w = cos x  cosh y  -  i sin x sinh y.
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    DEC       -10,+10      8400       4.5e-17     1.3e-17
+//    IEEE      -10,+10     30000       3.8e-16     1.0e-16
+
+// Cos returns the cosine of x.
+func Cos(x complex128) complex128 {
+	s, c := math.Sincos(real(x))
+	sh, ch := sinhcosh(imag(x))
+	return cmplx(c*ch, -s*sh)
+}
+
+// Complex hyperbolic cosine
+//
+// DESCRIPTION:
+//
+// ccosh(z) = cosh x  cos y + i sinh x sin y .
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    IEEE      -10,+10     30000       2.9e-16     8.1e-17
+
+// Cosh returns the hyperbolic cosine of x.
+func Cosh(x complex128) complex128 {
+	s, c := math.Sincos(imag(x))
+	sh, ch := sinhcosh(real(x))
+	return cmplx(c*ch, s*sh)
+}
+
+// calculate sinh and cosh
+func sinhcosh(x float64) (sh, ch float64) {
+	if math.Fabs(x) <= 0.5 {
+		return math.Sinh(x), math.Cosh(x)
+	}
+	e := math.Exp(x)
+	ei := 0.5 / e
+	e *= 0.5
+	return e - ei, e + ei
+}

src/pkg/cmath/sqrt.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// The original C code, the long comment, and the constants
+// below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
+// The go code is a simplified version of the original C.
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   moshier@na-net.ornl.gov
+
+// Complex square root
+//
+// DESCRIPTION:
+//
+// If z = x + iy,  r = |z|, then
+//
+//                       1/2
+// Re w  =  [ (r + x)/2 ]   ,
+//
+//                       1/2
+// Im w  =  [ (r - x)/2 ]   .
+//
+// Cancellation error in r-x or r+x is avoided by using the
+// identity  2 Re w Im w  =  y.
+//
+// Note that -w is also a square root of z.  The root chosen
+// is always in the right half plane and Im w has the same sign as y.
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    DEC       -10,+10     25000       3.2e-17     9.6e-18
+//    IEEE      -10,+10   1,000,000     2.9e-16     6.1e-17
+
+// Sqrt returns the square root of x.
+func Sqrt(x complex128) complex128 {
+	if imag(x) == 0 {
+		if real(x) == 0 {
+			return cmplx(0, 0)
+		}
+		if real(x) < 0 {
+			return cmplx(0, math.Sqrt(-real(x)))
+		}
+		return cmplx(math.Sqrt(real(x)), 0)
+	}
+	if real(x) == 0 {
+		if imag(x) < 0 {
+			r := math.Sqrt(-0.5 * imag(x))
+			return cmplx(r, -r)
+		}
+		r := math.Sqrt(0.5 * imag(x))
+		return cmplx(r, r)
+	}
+	a := real(x)
+	b := imag(x)
+	var scale float64
+	// Rescale to avoid internal overflow or underflow.
+	if math.Fabs(a) > 4 || math.Fabs(b) > 4 {
+		a *= 0.25
+		b *= 0.25
+		scale = 2
+	} else {
+		a *= 1.8014398509481984e16 // 2^54
+		b *= 1.8014398509481984e16
+		scale = 7.450580596923828125e-9 // 2^-27
+	}
+	r := math.Hypot(a, b)
+	var t float64
+	if a > 0 {
+		t = math.Sqrt(0.5*r + 0.5*a)
+		r = scale * math.Fabs((0.5*b)/t)
+		t *= scale
+	} else {
+		r = math.Sqrt(0.5*r - 0.5*a)
+		t = scale * math.Fabs((0.5*b)/r)
+		r *= scale
+	}
+	if b < 0 {
+		return cmplx(t, -r)
+	}
+	return cmplx(t, r)
+}

src/pkg/cmath/tan.go

+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cmath
+
+import "math"
+
+// The original C code, the long comment, and the constants
+// below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c.
+// The go code is a simplified version of the original C.
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   moshier@na-net.ornl.gov
+
+// Complex circular tangent
+//
+// DESCRIPTION:
+//
+// If
+//     z = x + iy,
+//
+// then
+//
+//           sin 2x  +  i sinh 2y
+//     w  =  --------------------.
+//            cos 2x  +  cosh 2y
+//
+// On the real axis the denominator is zero at odd multiples
+// of PI/2.  The denominator is evaluated by its Taylor
+// series near these points.
+//
+// ctan(z) = -i ctanh(iz).
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    DEC       -10,+10      5200       7.1e-17     1.6e-17
+//    IEEE      -10,+10     30000       7.2e-16     1.2e-16
+// Also tested by ctan * ccot = 1 and catan(ctan(z))  =  z.
+
+// Tan returns the tangent of x.
+func Tan(x complex128) complex128 {
+	d := math.Cos(2*real(x)) + math.Cosh(2*imag(x))
+	if math.Fabs(d) < 0.25 {
+		d = tanSeries(x)
+	}
+	if d == 0 {
+		return Inf()
+	}
+	return cmplx(math.Sin(2*real(x))/d, math.Sinh(2*imag(x))/d)
+}
+
+// Complex hyperbolic tangent
+//
+// DESCRIPTION:
+//
+// tanh z = (sinh 2x  +  i sin 2y) / (cosh 2x + cos 2y) .
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    IEEE      -10,+10     30000       1.7e-14     2.4e-16
+
+// Tanh returns the hyperbolic tangent of x.
+func Tanh(x complex128) complex128 {
+	d := math.Cosh(2*real(x)) + math.Cos(2*imag(x))
+	if d == 0 {
+		return Inf()
+	}
+	return cmplx(math.Sinh(2*real(x))/d, math.Sin(2*imag(x))/d)
+}
+
+// Program to subtract nearest integer multiple of PI
+func reducePi(x float64) float64 {
+	const (
+		// extended precision value of PI:
+		DP1 = 3.14159265160560607910E0   // ?? 0x400921fb54000000
+		DP2 = 1.98418714791870343106E-9  // ?? 0x3e210b4610000000
+		DP3 = 1.14423774522196636802E-17 // ?? 0x3c6a62633145c06e
+	)
+	t := x / math.Pi
+	if t >= 0 {
+		t += 0.5
+	} else {
+		t -= 0.5
+	}
+	t = float64(int64(t)) // int64(t) = the multiple
+	return ((x - t*DP1) - t*DP2) - t*DP3
+}
+
+// Taylor series expansion for cosh(2y) - cos(2x)
+func tanSeries(z complex128) float64 {
+	const MACHEP = 1.0 / (1 << 53)
+	x := math.Fabs(2 * real(z))
+	y := math.Fabs(2 * imag(z))
+	x = reducePi(x)
+	x = x * x
+	y = y * y
+	x2 := float64(1)
+	y2 := float64(1)
+	f := float64(1)
+	rn := float64(0)
+	d := float64(0)
+	for {
+		rn += 1
+		f *= rn
+		rn += 1
+		f *= rn
+		x2 *= x
+		y2 *= y
+		t := y2 + x2
+		t /= f
+		d += t
+
+		rn += 1
+		f *= rn
+		rn += 1
+		f *= rn
+		x2 *= x
+		y2 *= y
+		t = y2 - x2
+		t /= f
+		d += t
+		if math.Fabs(t/d) <= MACHEP {
+			break
+		}
+	}
+	return d
+}
+
+// Complex circular cotangent
+//
+// DESCRIPTION:
+//
+// If
+//     z = x + iy,
+//
+// then
+//
+//           sin 2x  -  i sinh 2y
+//     w  =  --------------------.
+//            cosh 2y  -  cos 2x
+//
+// On the real axis, the denominator has zeros at even
+// multiples of PI/2.  Near these points it is evaluated
+// by a Taylor series.
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    DEC       -10,+10      3000       6.5e-17     1.6e-17
+//    IEEE      -10,+10     30000       9.2e-16     1.2e-16
+// Also tested by ctan * ccot = 1 + i0.
+
+// Cot returns the cotangent of x.
+func Cot(x complex128) complex128 {
+	d := math.Cosh(2*imag(x)) - math.Cos(2*real(x))
+	if math.Fabs(d) < 0.25 {
+		d = tanSeries(x)
+	}
+	if d == 0 {
+		return Inf()
+	}
+	return cmplx(math.Sin(2*real(x))/d, -math.Sinh(2*imag(x))/d)
+}