1// Copyright 2009 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 math
6
7//extern sqrt
8func libc_sqrt(float64) float64
9
10func Sqrt(x float64) float64 {
11	return libc_sqrt(x)
12}
13
14// The original C code and the long comment below are
15// from FreeBSD's /usr/src/lib/msun/src/e_sqrt.c and
16// came with this notice. The go code is a simplified
17// version of the original C.
18//
19// ====================================================
20// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
21//
22// Developed at SunPro, a Sun Microsystems, Inc. business.
23// Permission to use, copy, modify, and distribute this
24// software is freely granted, provided that this notice
25// is preserved.
26// ====================================================
27//
28// __ieee754_sqrt(x)
29// Return correctly rounded sqrt.
30//           -----------------------------------------
31//           | Use the hardware sqrt if you have one |
32//           -----------------------------------------
33// Method:
34//   Bit by bit method using integer arithmetic. (Slow, but portable)
35//   1. Normalization
36//      Scale x to y in [1,4) with even powers of 2:
37//      find an integer k such that  1 <= (y=x*2**(2k)) < 4, then
38//              sqrt(x) = 2**k * sqrt(y)
39//   2. Bit by bit computation
40//      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
41//           i                                                   0
42//                                     i+1         2
43//          s  = 2*q , and      y  =  2   * ( y - q  ).          (1)
44//           i      i            i                 i
45//
46//      To compute q    from q , one checks whether
47//                  i+1       i
48//
49//                            -(i+1) 2
50//                      (q + 2      )  <= y.                     (2)
51//                        i
52//                                                            -(i+1)
53//      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
54//                             i+1   i             i+1   i
55//
56//      With some algebraic manipulation, it is not difficult to see
57//      that (2) is equivalent to
58//                             -(i+1)
59//                      s  +  2       <= y                       (3)
60//                       i                i
61//
62//      The advantage of (3) is that s  and y  can be computed by
63//                                    i      i
64//      the following recurrence formula:
65//          if (3) is false
66//
67//          s     =  s  ,       y    = y   ;                     (4)
68//           i+1      i          i+1    i
69//
70//      otherwise,
71//                         -i                      -(i+1)
72//          s     =  s  + 2  ,  y    = y  -  s  - 2              (5)
73//           i+1      i          i+1    i     i
74//
75//      One may easily use induction to prove (4) and (5).
76//      Note. Since the left hand side of (3) contain only i+2 bits,
77//            it is not necessary to do a full (53-bit) comparison
78//            in (3).
79//   3. Final rounding
80//      After generating the 53 bits result, we compute one more bit.
81//      Together with the remainder, we can decide whether the
82//      result is exact, bigger than 1/2ulp, or less than 1/2ulp
83//      (it will never equal to 1/2ulp).
84//      The rounding mode can be detected by checking whether
85//      huge + tiny is equal to huge, and whether huge - tiny is
86//      equal to huge for some floating point number "huge" and "tiny".
87//
88//
89// Notes:  Rounding mode detection omitted. The constants "mask", "shift",
90// and "bias" are found in src/math/bits.go
91
92// Sqrt returns the square root of x.
93//
94// Special cases are:
95//	Sqrt(+Inf) = +Inf
96//	Sqrt(±0) = ±0
97//	Sqrt(x < 0) = NaN
98//	Sqrt(NaN) = NaN
99
100// Note: Sqrt is implemented in assembly on some systems.
101// Others have assembly stubs that jump to func sqrt below.
102// On systems where Sqrt is a single instruction, the compiler
103// may turn a direct call into a direct use of that instruction instead.
104
105func sqrt(x float64) float64 {
106	// special cases
107	switch {
108	case x == 0 || IsNaN(x) || IsInf(x, 1):
109		return x
110	case x < 0:
111		return NaN()
112	}
113	ix := Float64bits(x)
114	// normalize x
115	exp := int((ix >> shift) & mask)
116	if exp == 0 { // subnormal x
117		for ix&(1<<shift) == 0 {
118			ix <<= 1
119			exp--
120		}
121		exp++
122	}
123	exp -= bias // unbias exponent
124	ix &^= mask << shift
125	ix |= 1 << shift
126	if exp&1 == 1 { // odd exp, double x to make it even
127		ix <<= 1
128	}
129	exp >>= 1 // exp = exp/2, exponent of square root
130	// generate sqrt(x) bit by bit
131	ix <<= 1
132	var q, s uint64               // q = sqrt(x)
133	r := uint64(1 << (shift + 1)) // r = moving bit from MSB to LSB
134	for r != 0 {
135		t := s + r
136		if t <= ix {
137			s = t + r
138			ix -= t
139			q += r
140		}
141		ix <<= 1
142		r >>= 1
143	}
144	// final rounding
145	if ix != 0 { // remainder, result not exact
146		q += q & 1 // round according to extra bit
147	}
148	ix = q>>1 + uint64(exp-1+bias)<<shift // significand + biased exponent
149	return Float64frombits(ix)
150}
151