Source file src/math/big/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 "math"
     8  
     9  var (
    10  	three = NewFloat(3.0)
    11  )
    12  
    13  // Sqrt sets z to the rounded square root of x, and returns it.
    14  //
    15  // If z's precision is 0, it is changed to x's precision before the
    16  // operation. Rounding is performed according to z's precision and
    17  // rounding mode.
    18  //
    19  // The function panics if z < 0. The value of z is undefined in that
    20  // case.
    21  func (z *Float) Sqrt(x *Float) *Float {
    22  	if debugFloat {
    23  		x.validate()
    24  	}
    25  
    26  	if z.prec == 0 {
    27  		z.prec = x.prec
    28  	}
    29  
    30  	if x.Sign() == -1 {
    31  		// following IEEE754-2008 (section 7.2)
    32  		panic(ErrNaN{"square root of negative operand"})
    33  	}
    34  
    35  	// handle ±0 and +∞
    36  	if x.form != finite {
    37  		z.acc = Exact
    38  		z.form = x.form
    39  		z.neg = x.neg // IEEE754-2008 requires √±0 = ±0
    40  		return z
    41  	}
    42  
    43  	// MantExp sets the argument's precision to the receiver's, and
    44  	// when z.prec > x.prec this will lower z.prec. Restore it after
    45  	// the MantExp call.
    46  	prec := z.prec
    47  	b := x.MantExp(z)
    48  	z.prec = prec
    49  
    50  	// Compute √(z·2**b) as
    51  	//   √( z)·2**(½b)     if b is even
    52  	//   √(2z)·2**(⌊½b⌋)   if b > 0 is odd
    53  	//   √(½z)·2**(⌈½b⌉)   if b < 0 is odd
    54  	switch b % 2 {
    55  	case 0:
    56  		// nothing to do
    57  	case 1:
    58  		z.exp++
    59  	case -1:
    60  		z.exp--
    61  	}
    62  	// 0.25 <= z < 2.0
    63  
    64  	// Solving x² - z = 0 directly requires a Quo call, but it's
    65  	// faster for small precisions.
    66  	//
    67  	// Solving 1/x² - z = 0 avoids the Quo call and is much faster for
    68  	// high precisions.
    69  	//
    70  	// 128bit precision is an empirically chosen threshold.
    71  	if z.prec <= 128 {
    72  		z.sqrtDirect(z)
    73  	} else {
    74  		z.sqrtInverse(z)
    75  	}
    76  
    77  	// re-attach halved exponent
    78  	return z.SetMantExp(z, b/2)
    79  }
    80  
    81  // Compute √x (up to prec 128) by solving
    82  //   t² - x = 0
    83  // for t, starting with a 53 bits precision guess from math.Sqrt and
    84  // then using at most two iterations of Newton's method.
    85  func (z *Float) sqrtDirect(x *Float) {
    86  	// let
    87  	//   f(t) = t² - x
    88  	// then
    89  	//   g(t) = f(t)/f'(t) = ½(t² - x)/t
    90  	// and the next guess is given by
    91  	//   t2 = t - g(t) = ½(t² + x)/t
    92  	u := new(Float)
    93  	ng := func(t *Float) *Float {
    94  		u.prec = t.prec
    95  		u.Mul(t, t)        // u = t²
    96  		u.Add(u, x)        //   = t² + x
    97  		u.exp--            //   = ½(t² + x)
    98  		return t.Quo(u, t) //   = ½(t² + x)/t
    99  	}
   100  
   101  	xf, _ := x.Float64()
   102  	sq := NewFloat(math.Sqrt(xf))
   103  
   104  	switch {
   105  	case z.prec > 128:
   106  		panic("sqrtDirect: only for z.prec <= 128")
   107  	case z.prec > 64:
   108  		sq.prec *= 2
   109  		sq = ng(sq)
   110  		fallthrough
   111  	default:
   112  		sq.prec *= 2
   113  		sq = ng(sq)
   114  	}
   115  
   116  	z.Set(sq)
   117  }
   118  
   119  // Compute √x (to z.prec precision) by solving
   120  //   1/t² - x = 0
   121  // for t (using Newton's method), and then inverting.
   122  func (z *Float) sqrtInverse(x *Float) {
   123  	// let
   124  	//   f(t) = 1/t² - x
   125  	// then
   126  	//   g(t) = f(t)/f'(t) = -½t(1 - xt²)
   127  	// and the next guess is given by
   128  	//   t2 = t - g(t) = ½t(3 - xt²)
   129  	u := newFloat(z.prec)
   130  	v := newFloat(z.prec)
   131  	ng := func(t *Float) *Float {
   132  		u.prec = t.prec
   133  		v.prec = t.prec
   134  		u.Mul(t, t)     // u = t²
   135  		u.Mul(x, u)     //   = xt²
   136  		v.Sub(three, u) // v = 3 - xt²
   137  		u.Mul(t, v)     // u = t(3 - xt²)
   138  		u.exp--         //   = ½t(3 - xt²)
   139  		return t.Set(u)
   140  
   141  	}
   142  
   143  	xf, _ := x.Float64()
   144  	sqi := newFloat(z.prec)
   145  	sqi.SetFloat64(1 / math.Sqrt(xf))
   146  	for prec := z.prec + 32; sqi.prec < prec; {
   147  		sqi.prec *= 2
   148  		sqi = ng(sqi)
   149  	}
   150  	// sqi = 1/√x
   151  
   152  	// x/√x = √x
   153  	z.Mul(x, sqi)
   154  }
   155  
   156  // newFloat returns a new *Float with space for twice the given
   157  // precision.
   158  func newFloat(prec2 uint32) *Float {
   159  	z := new(Float)
   160  	// nat.make ensures the slice length is > 0
   161  	z.mant = z.mant.make(int(prec2/_W) * 2)
   162  	return z
   163  }
   164  

View as plain text