Source file src/pkg/math/big/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 "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