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