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
		 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  //	 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