sqrt.go

     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