# Source file src/math/big/sqrt.go

```     1  // Copyright 2017 The Go Authors. All rights reserved.
2  // Use of this source code is governed by a BSD-style
4
5  package big
6
7  import (
8  	"math"
9  	"sync"
10  )
11
12  var threeOnce struct {
13  	sync.Once
14  	v *Float
15  }
16
17  func three() *Float {
18  	threeOnce.Do(func() {
19  		threeOnce.v = NewFloat(3.0)
20  	})
21  	return threeOnce.v
22  }
23
24  // Sqrt sets z to the rounded square root of x, and returns it.
25  //
26  // If z's precision is 0, it is changed to x's precision before the
27  // operation. Rounding is performed according to z's precision and
28  // rounding mode, but z's accuracy is not computed. Specifically, the
29  // result of z.Acc() is undefined.
30  //
31  // The function panics if z < 0. The value of z is undefined in that
32  // case.
33  func (z *Float) Sqrt(x *Float) *Float {
34  	if debugFloat {
35  		x.validate()
36  	}
37
38  	if z.prec == 0 {
39  		z.prec = x.prec
40  	}
41
42  	if x.Sign() == -1 {
43  		// following IEEE754-2008 (section 7.2)
44  		panic(ErrNaN{"square root of negative operand"})
45  	}
46
47  	// handle ±0 and +∞
48  	if x.form != finite {
49  		z.acc = Exact
50  		z.form = x.form
51  		z.neg = x.neg // IEEE754-2008 requires √±0 = ±0
52  		return z
53  	}
54
55  	// MantExp sets the argument's precision to the receiver's, and
56  	// when z.prec > x.prec this will lower z.prec. Restore it after
57  	// the MantExp call.
58  	prec := z.prec
59  	b := x.MantExp(z)
60  	z.prec = prec
61
62  	// Compute √(z·2**b) as
63  	//   √( z)·2**(½b)     if b is even
64  	//   √(2z)·2**(⌊½b⌋)   if b > 0 is odd
65  	//   √(½z)·2**(⌈½b⌉)   if b < 0 is odd
66  	switch b % 2 {
67  	case 0:
68  		// nothing to do
69  	case 1:
70  		z.exp++
71  	case -1:
72  		z.exp--
73  	}
74  	// 0.25 <= z < 2.0
75
76  	// Solving 1/x² - z = 0 avoids Quo calls and is faster, especially
77  	// for high precisions.
78  	z.sqrtInverse(z)
79
80  	// re-attach halved exponent
81  	return z.SetMantExp(z, b/2)
82  }
83
84  // Compute √x (to z.prec precision) by solving
85  //
86  //	1/t² - x = 0
87  //
88  // for t (using Newton's method), and then inverting.
89  func (z *Float) sqrtInverse(x *Float) {
90  	// let
91  	//   f(t) = 1/t² - x
92  	// then
93  	//   g(t) = f(t)/f'(t) = -½t(1 - xt²)
94  	// and the next guess is given by
95  	//   t2 = t - g(t) = ½t(3 - xt²)
96  	u := newFloat(z.prec)
97  	v := newFloat(z.prec)
98  	three := three()
99  	ng := func(t *Float) *Float {
100  		u.prec = t.prec
101  		v.prec = t.prec
102  		u.Mul(t, t)     // u = t²
103  		u.Mul(x, u)     //   = xt²
104  		v.Sub(three, u) // v = 3 - xt²
105  		u.Mul(t, v)     // u = t(3 - xt²)
106  		u.exp--         //   = ½t(3 - xt²)
107  		return t.Set(u)
108  	}
109
110  	xf, _ := x.Float64()
111  	sqi := newFloat(z.prec)
112  	sqi.SetFloat64(1 / math.Sqrt(xf))
113  	for prec := z.prec + 32; sqi.prec < prec; {
114  		sqi.prec *= 2
115  		sqi = ng(sqi)
116  	}
117  	// sqi = 1/√x
118
119  	// x/√x = √x
120  	z.Mul(x, sqi)
121  }
122
123  // newFloat returns a new *Float with space for twice the given
124  // precision.
125  func newFloat(prec2 uint32) *Float {
126  	z := new(Float)
127  	// nat.make ensures the slice length is > 0
128  	z.mant = z.mant.make(int(prec2/_W) * 2)
129  	return z
130  }
131
```

View as plain text