...
Run Format

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

View as plain text