````// Copyright 2009 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 math`

`// The original C code and the long comment below are`
`// from FreeBSD's /usr/src/lib/msun/src/e_sqrt.c and`
`// came with this notice. The go code is a simplified`
`// version of the original C.`
`//`
`// ====================================================`
`// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.`
`//`
`// Developed at SunPro, a Sun Microsystems, Inc. business.`
`// Permission to use, copy, modify, and distribute this`
`// software is freely granted, provided that this notice`
`// is preserved.`
`// ====================================================`
`//`
`// __ieee754_sqrt(x)`
`// Return correctly rounded sqrt.`
`//           -----------------------------------------`
`//           | Use the hardware sqrt if you have one |`
`//           -----------------------------------------`
`// Method:`
`//   Bit by bit method using integer arithmetic. (Slow, but portable)`
`//   1. Normalization`
`//      Scale x to y in [1,4) with even powers of 2:`
`//      find an integer k such that  1 <= (y=x*2**(2k)) < 4, then`
`//              sqrt(x) = 2**k * sqrt(y)`
`//   2. Bit by bit computation`
`//      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),`
`//           i                                                   0`
`//                                     i+1         2`
`//          s  = 2*q , and      y  =  2   * ( y - q  ).          (1)`
`//           i      i            i                 i`
`//`
`//      To compute q    from q , one checks whether`
`//                  i+1       i`
`//`
`//                            -(i+1) 2`
`//                      (q + 2      )  <= y.                     (2)`
`//                        i`
`//                                                            -(i+1)`
`//      If (2) is false, then q   = q ; otherwise q   = q  + 2      .`
`//                             i+1   i             i+1   i`
`//`
`//      With some algebraic manipulation, it is not difficult to see`
`//      that (2) is equivalent to`
`//                             -(i+1)`
`//                      s  +  2       <= y                       (3)`
`//                       i                i`
`//`
`//      The advantage of (3) is that s  and y  can be computed by`
`//                                    i      i`
`//      the following recurrence formula:`
`//          if (3) is false`
`//`
`//          s     =  s  ,       y    = y   ;                     (4)`
`//           i+1      i          i+1    i`
`//`
`//      otherwise,`
`//                         -i                      -(i+1)`
`//          s     =  s  + 2  ,  y    = y  -  s  - 2              (5)`
`//           i+1      i          i+1    i     i`
`//`
`//      One may easily use induction to prove (4) and (5).`
`//      Note. Since the left hand side of (3) contain only i+2 bits,`
`//            it is not necessary to do a full (53-bit) comparison`
`//            in (3).`
`//   3. Final rounding`
`//      After generating the 53 bits result, we compute one more bit.`
`//      Together with the remainder, we can decide whether the`
`//      result is exact, bigger than 1/2ulp, or less than 1/2ulp`
`//      (it will never equal to 1/2ulp).`
`//      The rounding mode can be detected by checking whether`
`//      huge + tiny is equal to huge, and whether huge - tiny is`
`//      equal to huge for some floating point number "huge" and "tiny".`
`//`
`//`
`// Notes:  Rounding mode detection omitted. The constants "mask", "shift",`
`// and "bias" are found in src/math/bits.go`

`// Sqrt returns the square root of x.`
`//`
`// Special cases are:`
`//	Sqrt(+Inf) = +Inf`
`//	Sqrt(±0) = ±0`
`//	Sqrt(x < 0) = NaN`
`//	Sqrt(NaN) = NaN`
`func Sqrt(x float64) float64 {`
`	if haveArchSqrt {`
`		return archSqrt(x)`
`	}`
`	return sqrt(x)`
`}`

`// Note: Sqrt is implemented in assembly on some systems.`
`// Others have assembly stubs that jump to func sqrt below.`
`// On systems where Sqrt is a single instruction, the compiler`
`// may turn a direct call into a direct use of that instruction instead.`

`func sqrt(x float64) float64 {`
`	// special cases`
`	switch {`
`	case x == 0 || IsNaN(x) || IsInf(x, 1):`
`		return x`
`	case x < 0:`
`		return NaN()`
`	}`
`	ix := Float64bits(x)`
`	// normalize x`
`	exp := int((ix >> shift) & mask)`
`	if exp == 0 { // subnormal x`
`		for ix&(1<<shift) == 0 {`
`			ix <<= 1`
`			exp--`
`		}`
`		exp++`
`	}`
`	exp -= bias // unbias exponent`
`	ix &^= mask << shift`
`	ix |= 1 << shift`
`	if exp&1 == 1 { // odd exp, double x to make it even`
`		ix <<= 1`
`	}`
`	exp >>= 1 // exp = exp/2, exponent of square root`
`	// generate sqrt(x) bit by bit`
`	ix <<= 1`
`	var q, s uint64               // q = sqrt(x)`
`	r := uint64(1 << (shift + 1)) // r = moving bit from MSB to LSB`
`	for r != 0 {`
`		t := s + r`
`		if t <= ix {`
`			s = t + r`
`			ix -= t`
`			q += r`
`		}`
`		ix <<= 1`
`		r >>= 1`
`	}`
`	// final rounding`
`	if ix != 0 { // remainder, result not exact`
`		q += q & 1 // round according to extra bit`
`	}`
`	ix = q>>1 + uint64(exp-1+bias)<<shift // significand + biased exponent`
`	return Float64frombits(ix)`
`}`
```