...

Source file src/math/cmplx/sqrt.go

Documentation: math/cmplx

		 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  //	 [email protected]
		30  
		31  // Complex square root
		32  //
		33  // DESCRIPTION:
		34  //
		35  // If z = x + iy,	r = |z|, then
		36  //
		37  //											 1/2
		38  // Re w	=	[ (r + x)/2 ]	 ,
		39  //
		40  //											 1/2
		41  // Im w	=	[ (r - x)/2 ]	 .
		42  //
		43  // Cancellation error in r-x or r+x is avoided by using the
		44  // identity	2 Re w Im w	=	y.
		45  //
		46  // Note that -w is also a square root of z. The root chosen
		47  // is always in the right half plane and Im w has the same sign as y.
		48  //
		49  // ACCURACY:
		50  //
		51  //											Relative error:
		52  // arithmetic	 domain		 # trials			peak				 rms
		53  //		DEC			 -10,+10		 25000			 3.2e-17		 9.6e-18
		54  //		IEEE			-10,+10	 1,000,000		 2.9e-16		 6.1e-17
		55  
		56  // Sqrt returns the square root of x.
		57  // The result r is chosen so that real(r) ≥ 0 and imag(r) has the same sign as imag(x).
		58  func Sqrt(x complex128) complex128 {
		59  	if imag(x) == 0 {
		60  		// Ensure that imag(r) has the same sign as imag(x) for imag(x) == signed zero.
		61  		if real(x) == 0 {
		62  			return complex(0, imag(x))
		63  		}
		64  		if real(x) < 0 {
		65  			return complex(0, math.Copysign(math.Sqrt(-real(x)), imag(x)))
		66  		}
		67  		return complex(math.Sqrt(real(x)), imag(x))
		68  	} else if math.IsInf(imag(x), 0) {
		69  		return complex(math.Inf(1.0), imag(x))
		70  	}
		71  	if real(x) == 0 {
		72  		if imag(x) < 0 {
		73  			r := math.Sqrt(-0.5 * imag(x))
		74  			return complex(r, -r)
		75  		}
		76  		r := math.Sqrt(0.5 * imag(x))
		77  		return complex(r, r)
		78  	}
		79  	a := real(x)
		80  	b := imag(x)
		81  	var scale float64
		82  	// Rescale to avoid internal overflow or underflow.
		83  	if math.Abs(a) > 4 || math.Abs(b) > 4 {
		84  		a *= 0.25
		85  		b *= 0.25
		86  		scale = 2
		87  	} else {
		88  		a *= 1.8014398509481984e16 // 2**54
		89  		b *= 1.8014398509481984e16
		90  		scale = 7.450580596923828125e-9 // 2**-27
		91  	}
		92  	r := math.Hypot(a, b)
		93  	var t float64
		94  	if a > 0 {
		95  		t = math.Sqrt(0.5*r + 0.5*a)
		96  		r = scale * math.Abs((0.5*b)/t)
		97  		t *= scale
		98  	} else {
		99  		r = math.Sqrt(0.5*r - 0.5*a)
	 100  		t = scale * math.Abs((0.5*b)/r)
	 101  		r *= scale
	 102  	}
	 103  	if b < 0 {
	 104  		return complex(t, -r)
	 105  	}
	 106  	return complex(t, r)
	 107  }
	 108  

View as plain text