codekingpro/portable-devtools

1// Copyright 2017 The Go Authors. All rights reserved.2// Use of this source code is governed by a BSD-style3// license that can be found in the LICENSE file.4 5package big6 7import (8	"math"9	"sync"10)11 12var threeOnce struct {13	sync.Once14	v *Float15}16 17func three() *Float {18	threeOnce.Do(func() {19		threeOnce.v = NewFloat(3.0)20	})21	return threeOnce.v22}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 the27// operation. Rounding is performed according to z's precision and28// rounding mode, but z's accuracy is not computed. Specifically, the29// result of z.Acc() is undefined.30//31// The function panics if z < 0. The value of z is undefined in that32// 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.prec40	}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 = Exact50		z.form = x.form51		z.neg = x.neg // IEEE754-2008 requires √±0 = ±052		return z53	}54 55	// MantExp sets the argument's precision to the receiver's, and56	// when z.prec > x.prec this will lower z.prec. Restore it after57	// the MantExp call.58	prec := z.prec59	b := x.MantExp(z)60	z.prec = prec61 62	// Compute √(z·2**b) as63	//   √( z)·2**(½b)     if b is even64	//   √(2z)·2**(⌊½b⌋)   if b > 0 is odd65	//   √(½z)·2**(⌈½b⌉)   if b < 0 is odd66	switch b % 2 {67	case 0:68		// nothing to do69	case 1:70		z.exp++71	case -1:72		z.exp--73	}74	// 0.25 <= z < 2.075 76	// Solving 1/x² - z = 0 avoids Quo calls and is faster, especially77	// for high precisions.78	z.sqrtInverse(z)79 80	// re-attach halved exponent81	return z.SetMantExp(z, b/2)82}83 84// Compute √x (to z.prec precision) by solving85//86//	1/t² - x = 087//88// for t (using Newton's method), and then inverting.89func (z *Float) sqrtInverse(x *Float) {90	// let91	//   f(t) = 1/t² - x92	// then93	//   g(t) = f(t)/f'(t) = -½t(1 - xt²)94	// and the next guess is given by95	//   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.prec101		v.prec = t.prec102		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 *= 2115		sqi = ng(sqi)116	}117	// sqi = 1/√x118 119	// x/√x = √x120	z.Mul(x, sqi)121}122 123// newFloat returns a new *Float with space for twice the given124// precision.125func newFloat(prec2 uint32) *Float {126	z := new(Float)127	// nat.make ensures the slice length is > 0128	z.mant = z.mant.make(int(prec2/_W) * 2)129	return z130}131