sqrt.go: remove allocations from Newton functions

name               old time/op    new time/op    delta
Exp2Prec53-4         69.4µs ± 4%    61.0µs ± 1%  -12.09%          (p=0.008 n=5+5)
Exp2Prec100-4         154µs ± 0%     137µs ± 0%  -11.39%          (p=0.008 n=5+5)
Exp2Prec1000-4       1.07ms ± 0%    0.98ms ± 0%   -7.91%          (p=0.008 n=5+5)
Exp2Prec10000-4      14.5ms ± 0%    14.3ms ± 0%   -0.97%          (p=0.008 n=5+5)
Exp2Prec100000-4      426ms ± 0%     426ms ± 0%     ~             (p=0.222 n=5+5)
Log2Prec53-4         63.1µs ± 0%    56.4µs ± 0%  -10.63%          (p=0.008 n=5+5)
Log2Prec100-4        69.9µs ± 0%    64.7µs ± 3%   -7.44%          (p=0.008 n=5+5)
Log2Prec1000-4        216µs ± 0%     196µs ± 1%   -9.02%          (p=0.008 n=5+5)
Log2Prec10000-4      3.65ms ± 0%    3.60ms ± 0%   -1.33%          (p=0.008 n=5+5)
Log2Prec100000-4      134ms ± 2%     133ms ± 0%     ~             (p=0.151 n=5+5)
Sqrt2Prec53-4        1.16µs ± 0%    1.16µs ± 0%     ~             (p=0.087 n=5+5)
Sqrt2Prec100-4       2.38µs ± 4%    2.27µs ± 5%     ~             (p=0.127 n=5+5)
Sqrt2Prec1000-4      8.47µs ± 0%    7.28µs ± 0%  -14.07%          (p=0.008 n=5+5)
Sqrt2Prec10000-4     51.4µs ± 1%    49.6µs ± 0%   -3.43%          (p=0.008 n=5+5)
Sqrt2Prec100000-4     972µs ± 0%     974µs ± 0%     ~             (p=0.151 n=5+5)

name               old allocs/op  new allocs/op  delta
Exp2Prec53-4            593 ± 0%       461 ± 0%  -22.26%          (p=0.008 n=5+5)
Exp2Prec100-4         1.28k ± 0%     0.99k ± 0%  -22.46%          (p=0.008 n=5+5)
Exp2Prec1000-4        5.86k ± 0%     4.35k ± 0%  -25.76%          (p=0.008 n=5+5)
Exp2Prec10000-4       13.5k ± 0%      9.9k ± 0%  -26.90%          (p=0.008 n=5+5)
Exp2Prec100000-4      27.0k ± 0%     19.7k ± 0%  -27.19%          (p=0.008 n=5+5)
Log2Prec53-4            539 ± 0%       419 ± 0%  -22.26%          (p=0.008 n=5+5)
Log2Prec100-4           590 ± 0%       458 ± 0%  -22.37%          (p=0.008 n=5+5)
Log2Prec1000-4        1.24k ± 0%     0.90k ± 0%  -27.53%          (p=0.008 n=5+5)
Log2Prec10000-4       2.70k ± 0%     1.96k ± 0%  -27.28%          (p=0.008 n=5+5)
Log2Prec100000-4      5.21k ± 0%     3.89k ± 0%  -25.35%          (p=0.008 n=5+5)
Sqrt2Prec53-4          14.0 ± 0%      14.0 ± 0%     ~     (all samples are equal)
Sqrt2Prec100-4         23.0 ± 0%      21.0 ± 0%   -8.70%          (p=0.008 n=5+5)
Sqrt2Prec1000-4        63.0 ± 0%      43.0 ± 0%  -31.75%          (p=0.008 n=5+5)
Sqrt2Prec10000-4       92.0 ± 0%      60.0 ± 0%  -34.78%          (p=0.008 n=5+5)
Sqrt2Prec100000-4       128 ± 0%        84 ± 0%  -34.38%          (p=0.008 n=5+5)

sqrt.go

@@ -1,15 +1,16 @@
-// Package floats provides the implementation of a few additional operations for the
-// standard library big.Float type.
+// Package floats provides the implementation of a few additional
+// operations for the standard library big.Float type.
 package floats
 
 import (
 	"math"
 	"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