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