Black Lives Matter. Support the Equal Justice Initiative.

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

## Documentation: math/big

```     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  //   1/t² - x = 0
86  // for t (using Newton's method), and then inverting.
87  func (z *Float) sqrtInverse(x *Float) {
88  	// let
89  	//   f(t) = 1/t² - x
90  	// then
91  	//   g(t) = f(t)/f'(t) = -½t(1 - xt²)
92  	// and the next guess is given by
93  	//   t2 = t - g(t) = ½t(3 - xt²)
94  	u := newFloat(z.prec)
95  	v := newFloat(z.prec)
96  	three := three()
97  	ng := func(t *Float) *Float {
98  		u.prec = t.prec
99  		v.prec = t.prec
100  		u.Mul(t, t)     // u = t²
101  		u.Mul(x, u)     //   = xt²
102  		v.Sub(three, u) // v = 3 - xt²
103  		u.Mul(t, v)     // u = t(3 - xt²)
104  		u.exp--         //   = ½t(3 - xt²)
105  		return t.Set(u)
106  	}
107
108  	xf, _ := x.Float64()
109  	sqi := newFloat(z.prec)
110  	sqi.SetFloat64(1 / math.Sqrt(xf))
111  	for prec := z.prec + 32; sqi.prec < prec; {
112  		sqi.prec *= 2
113  		sqi = ng(sqi)
114  	}
115  	// sqi = 1/√x
116
117  	// x/√x = √x
118  	z.Mul(x, sqi)
119  }
120
121  // newFloat returns a new *Float with space for twice the given
122  // precision.
123  func newFloat(prec2 uint32) *Float {
124  	z := new(Float)
125  	// nat.make ensures the slice length is > 0
126  	z.mant = z.mant.make(int(prec2/_W) * 2)
127  	return z
128  }
129
```

View as plain text