Merge branch 'dev'

sqrt.go

@@ -7,9 +7,10 @@ import (
 	"math/big"
 )
 
-// Sqrt returns a big.Float representation of the square root of z. Precision is
-// the same as the one of the argument. The function panics if z is negative, returns ±0
-// when z = ±0, and +Inf when z = +Inf.
+// Sqrt returns a big.Float representation of the square root of
+// z. Precision is the same as the one of the argument. The function
+// panics if z is negative, returns ±0 when z = ±0, and +Inf when z =
+// +Inf.
 func Sqrt(z *big.Float) *big.Float {
 
 	// panic on negative z
@@ -32,9 +33,8 @@ func Sqrt(z *big.Float) *big.Float {
 	//   √(2a)·2**b/2      if b > 0 is odd
 	//   √(0.5a)·2**b/2    if b < 0 is odd
 	//
-	// The difference in the odd exponent case is due
-	// to the fact that exp/2 is rounded in different
-	// directions when exp is negative.
+	// The difference in the odd exponent case is due to the fact that
+	// exp/2 is rounded in different directions when exp is negative.
 	mant := new(big.Float)
 	exp := z.MantExp(mant)
 	switch exp % 2 {
@@ -44,12 +44,13 @@ func Sqrt(z *big.Float) *big.Float {
 		mant.Mul(big.NewFloat(0.5), mant)
 	}
 
-	// Solving x² - z = 0 directly requires a Quo
-	// call, but it's faster for small precisions.
-	// Solvin 1/x² - z = 0 avoids the Quo call and
-	// is much faster for high precisions.
-	// Use sqrtDirect for prec <= 128 and
-	// sqrtInverse for prec > 128.
+	// Solving x² - z = 0 directly requires a Quo call, but it's
+	// faster for small precisions.
+	//
+	// Solving 1/x² - z = 0 avoids the Quo call and is much faster for
+	// high precisions.
+	//
+	// Use sqrtDirect for prec <= 128 and sqrtInverse for prec > 128.
 	var x *big.Float
 	if z.Prec() <= 128 {
 		x = sqrtDirect(mant)
@@ -66,10 +67,11 @@ func Sqrt(z *big.Float) *big.Float {
 // t² - z = 0 for t
 func sqrtDirect(z *big.Float) *big.Float {
 	// f(t)/f'(t) = 0.5(t² - z)/t
+	half := big.NewFloat(0.5)
 	f := func(t *big.Float) *big.Float {
 		x := new(big.Float).Mul(t, t) // x = t²
 		x.Sub(x, z)                   // x = t² - z
-		x.Mul(big.NewFloat(0.5), x)   // x = 0.5(t² - z)
+		x.Mul(half, x)                // x = 0.5(t² - z)
 		return x.Quo(x, t)            // return x = 0.5(t² - z)/t
 	}
 
@@ -84,13 +86,15 @@ func sqrtDirect(z *big.Float) *big.Float {
 // 1/t² - z = 0 for x and then inverting.
 func sqrtInverse(z *big.Float) *big.Float {
 	// f(t)/f'(t) = -0.5t(1 - zt²)
+	nhalf := big.NewFloat(-0.5)
+	one := big.NewFloat(1)
 	f := func(t *big.Float) *big.Float {
 		u := new(big.Float)
 		u.Mul(t, t)                     // u = t²
 		u.Mul(u, z)                     // u = zt²
-		u.Sub(big.NewFloat(1), u)       // u = 1 - zt²
-		u.Mul(u, big.NewFloat(-0.5))    // u = 0.5(1 - zt²)
-		return new(big.Float).Mul(t, u) // x = 0.5t(1 - zt²)
+		u.Sub(one, u)                   // u = 1 - zt²
+		u.Mul(u, nhalf)                 // u = -0.5(1 - zt²)
+		return new(big.Float).Mul(t, u) // x = -0.5t(1 - zt²)
 	}
 
 	// initial guess