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 (
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 := &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