Source file src/pkg/math/sqrt.go

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