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
