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
5package big
6
7import (
8	"math"
9	"sync"
10)
11
12var threeOnce struct {
13	sync.Once
14	v *Float
15}
16
17func 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.
33func (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.
89func (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.
125func 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