Home | History | Annotate | Download | only in runtime
      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