From 856dccb1757b71d4d1a85270db79d42120cdf0e9 Mon Sep 17 00:00:00 2001
From: Alberto Donizetti <alb.donizetti@gmail.com>
Date: Wed, 1 Nov 2017 12:36:56 +0100
Subject: [PATCH] math/big: avoid unnecessary Newton iteration in Float.Sqrt
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

An initial draft of the Newton code for Float.Sqrt was structured like
this:

  for condition
    // do Newton iteration..
    prec *= 2

since prec, at the end of the loop, was double the precision used in
the last Newton iteration, the termination condition was set to
2*limit. The code was later rewritten in the form

  for condition
    prec *= 2
    // do Newton iteration..

but condition was not updated, and it's still 2*limit, which is about
double what we actually need, and is triggering the execution of an
additional, and unnecessary, Newton iteration.

This change adjusts the Newton termination condition to the (correct)
value of z.prec, plus 32 guard bits as a safety margin.

name                 old time/op    new time/op    delta
FloatSqrt/64-4          798ns 卤 3%     802ns 卤 3%     ~     (p=0.458 n=8+8)
FloatSqrt/128-4        1.65碌s 卤 1%    1.65碌s 卤 1%     ~     (p=0.290 n=8+8)
FloatSqrt/256-4        3.10碌s 卤 1%    2.10碌s 卤 0%  -32.32%  (p=0.000 n=8+7)
FloatSqrt/1000-4       8.83碌s 卤 1%    4.91碌s 卤 2%  -44.39%  (p=0.000 n=8+8)
FloatSqrt/10000-4       107碌s 卤 1%      40碌s 卤 1%  -62.68%  (p=0.000 n=8+8)
FloatSqrt/100000-4     2.91ms 卤 1%    0.96ms 卤 1%  -67.13%  (p=0.000 n=8+8)
FloatSqrt/1000000-4     240ms 卤 1%      80ms 卤 1%  -66.66%  (p=0.000 n=8+8)

name                 old alloc/op   new alloc/op   delta
FloatSqrt/64-4           416B 卤 0%      416B 卤 0%     ~     (all equal)
FloatSqrt/128-4          720B 卤 0%      720B 卤 0%     ~     (all equal)
FloatSqrt/256-4        1.34kB 卤 0%    0.82kB 卤 0%  -39.29%  (p=0.000 n=8+8)
FloatSqrt/1000-4       5.09kB 卤 0%    2.50kB 卤 0%  -50.94%  (p=0.000 n=8+8)
FloatSqrt/10000-4      45.9kB 卤 0%    23.5kB 卤 0%  -48.81%  (p=0.000 n=8+8)
FloatSqrt/100000-4      533kB 卤 0%     251kB 卤 0%  -52.90%  (p=0.000 n=8+8)
FloatSqrt/1000000-4    9.21MB 卤 0%    4.61MB 卤 0%  -49.98%  (p=0.000 n=8+8)

name                 old allocs/op  new allocs/op  delta
FloatSqrt/64-4           9.00 卤 0%      9.00 卤 0%     ~     (all equal)
FloatSqrt/128-4          13.0 卤 0%      13.0 卤 0%     ~     (all equal)
FloatSqrt/256-4          15.0 卤 0%      12.0 卤 0%  -20.00%  (p=0.000 n=8+8)
FloatSqrt/1000-4         24.0 卤 0%      19.0 卤 0%  -20.83%  (p=0.000 n=8+8)
FloatSqrt/10000-4        40.0 卤 0%      35.0 卤 0%  -12.50%  (p=0.000 n=8+8)
FloatSqrt/100000-4       66.0 卤 0%      55.0 卤 0%  -16.67%  (p=0.000 n=8+8)
FloatSqrt/1000000-4       143 卤 0%       122 卤 0%  -14.69%  (p=0.000 n=8+8)

Change-Id: I4868adb7f8960f2ca20e7792734c2e6211669fc0
Reviewed-on: https://go-review.googlesource.com/75010
Reviewed-by: Robert Griesemer <gri@golang.org>
---
 src/math/big/sqrt.go | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/math/big/sqrt.go b/src/math/big/sqrt.go
index d1deb77652..00433cfe7a 100644
--- a/src/math/big/sqrt.go
+++ b/src/math/big/sqrt.go
@@ -140,7 +140,7 @@ func (z *Float) sqrtInverse(x *Float) {
 
 	xf, _ := x.Float64()
 	sqi := NewFloat(1 / math.Sqrt(xf))
-	for prec := 2 * z.prec; sqi.prec < prec; {
+	for prec := z.prec + 32; sqi.prec < prec; {
 		sqi.prec *= 2
 		sqi = ng(sqi)
 	}
-- 
2.30.9