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