1 // Copyright 2017 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 package big 6 7 import "math" 8 9 var ( 10 half = NewFloat(0.5) 11 two = NewFloat(2.0) 12 three = NewFloat(3.0) 13 ) 14 15 // Sqrt sets z to the rounded square root of x, and returns it. 16 // 17 // If z's precision is 0, it is changed to x's precision before the 18 // operation. Rounding is performed according to z's precision and 19 // rounding mode. 20 // 21 // The function panics if z < 0. The value of z is undefined in that 22 // case. 23 func (z *Float) Sqrt(x *Float) *Float { 24 if debugFloat { 25 x.validate() 26 } 27 28 if z.prec == 0 { 29 z.prec = x.prec 30 } 31 32 if x.Sign() == -1 { 33 // following IEEE754-2008 (section 7.2) 34 panic(ErrNaN{"square root of negative operand"}) 35 } 36 37 // handle 0 and + 38 if x.form != finite { 39 z.acc = Exact 40 z.form = x.form 41 z.neg = x.neg // IEEE754-2008 requires 0 = 0 42 return z 43 } 44 45 // MantExp sets the argument's precision to the receiver's, and 46 // when z.prec > x.prec this will lower z.prec. Restore it after 47 // the MantExp call. 48 prec := z.prec 49 b := x.MantExp(z) 50 z.prec = prec 51 52 // Compute (z2**b) as 53 // ( z)2**(b) if b is even 54 // (2z)2**(b) if b > 0 is odd 55 // (z)2**(b) if b < 0 is odd 56 switch b % 2 { 57 case 0: 58 // nothing to do 59 case 1: 60 z.Mul(two, z) 61 case -1: 62 z.Mul(half, z) 63 } 64 // 0.25 <= z < 2.0 65 66 // Solving x - z = 0 directly requires a Quo call, but it's 67 // faster for small precisions. 68 // 69 // Solving 1/x - z = 0 avoids the Quo call and is much faster for 70 // high precisions. 71 // 72 // 128bit precision is an empirically chosen threshold. 73 if z.prec <= 128 { 74 z.sqrtDirect(z) 75 } else { 76 z.sqrtInverse(z) 77 } 78 79 // re-attach halved exponent 80 return z.SetMantExp(z, b/2) 81 } 82 83 // Compute x (up to prec 128) by solving 84 // t - x = 0 85 // for t, starting with a 53 bits precision guess from math.Sqrt and 86 // then using at most two iterations of Newton's method. 87 func (z *Float) sqrtDirect(x *Float) { 88 // let 89 // f(t) = t - x 90 // then 91 // g(t) = f(t)/f'(t) = (t - x)/t 92 // and the next guess is given by 93 // t2 = t - g(t) = (t + x)/t 94 u := new(Float) 95 ng := func(t *Float) *Float { 96 u.prec = t.prec 97 u.Mul(t, t) // u = t 98 u.Add(u, x) // = t + x 99 u.Mul(half, u) // = (t + x) 100 return t.Quo(u, t) // = (t + x)/t 101 } 102 103 xf, _ := x.Float64() 104 sq := NewFloat(math.Sqrt(xf)) 105 106 switch { 107 case z.prec > 128: 108 panic("sqrtDirect: only for z.prec <= 128") 109 case z.prec > 64: 110 sq.prec *= 2 111 sq = ng(sq) 112 fallthrough 113 default: 114 sq.prec *= 2 115 sq = ng(sq) 116 } 117 118 z.Set(sq) 119 } 120 121 // Compute x (to z.prec precision) by solving 122 // 1/t - x = 0 123 // for t (using Newton's method), and then inverting. 124 func (z *Float) sqrtInverse(x *Float) { 125 // let 126 // f(t) = 1/t - x 127 // then 128 // g(t) = f(t)/f'(t) = -t(1 - xt) 129 // and the next guess is given by 130 // t2 = t - g(t) = t(3 - xt) 131 u := new(Float) 132 ng := func(t *Float) *Float { 133 u.prec = t.prec 134 u.Mul(t, t) // u = t 135 u.Mul(x, u) // = xt 136 u.Sub(three, u) // = 3 - xt 137 u.Mul(t, u) // = t(3 - xt) 138 return t.Mul(half, u) // = t(3 - xt) 139 } 140 141 xf, _ := x.Float64() 142 sqi := NewFloat(1 / math.Sqrt(xf)) 143 for prec := z.prec + 32; sqi.prec < prec; { 144 sqi.prec *= 2 145 sqi = ng(sqi) 146 } 147 // sqi = 1/x 148 149 // x/x = x 150 z.Mul(x, sqi) 151 } 152
