...
Run Format

Source file src/runtime/sqrt.go

  // 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.
  
  // Copy of math/sqrt.go, here for use by ARM softfloat.
  // Modified to not use any floating point arithmetic so
  // that we don't clobber any floating-point registers
  // while emulating the sqrt instruction.
  
  package runtime
  
  // 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 does 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.
  
  const (
  	float64Mask  = 0x7FF
  	float64Shift = 64 - 11 - 1
  	float64Bias  = 1023
  	float64NaN   = 0x7FF8000000000001
  	float64Inf   = 0x7FF0000000000000
  	maxFloat64   = 1.797693134862315708145274237317043567981e+308 // 2**1023 * (2**53 - 1) / 2**52
  )
  
  // isnanu returns whether ix represents a NaN floating point number.
  func isnanu(ix uint64) bool {
  	exp := (ix >> float64Shift) & float64Mask
  	sig := ix << (64 - float64Shift) >> (64 - float64Shift)
  	return exp == float64Mask && sig != 0
  }
  
  func sqrt(ix uint64) uint64 {
  	// special cases
  	switch {
  	case ix == 0 || ix == 1<<63: // x == 0
  		return ix
  	case isnanu(ix): // x != x
  		return ix
  	case ix&(1<<63) != 0: // x < 0
  		return float64NaN
  	case ix == float64Inf: // x > MaxFloat
  		return ix
  	}
  	// normalize x
  	exp := int((ix >> float64Shift) & float64Mask)
  	if exp == 0 { // subnormal x
  		for ix&(1<<float64Shift) == 0 {
  			ix <<= 1
  			exp--
  		}
  		exp++
  	}
  	exp -= float64Bias // unbias exponent
  	ix &^= float64Mask << float64Shift
  	ix |= 1 << float64Shift
  	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 << (float64Shift + 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+float64Bias)<<float64Shift // significand + biased exponent
  	return ix
  }
  

View as plain text