平方根算法

笔者最近在看《计算机程序的构造和解释》一书，书中第一章讲到了平方根计算算法，笔者当时就在想一些脑中的平方根算法，就写了本文。

如果不谈论数学，工程层面上，求取一个平方根的实质是在限定的潜在解空间内搜索一个符合要求的值，潜在的值按照大小排列。最简单直白的就是使用二分的策略：假设要求数X的平方根，实质上可以化简为求数abs(X-Y*Y) <= N(N为常数)，Y的取值范围[0, X]。这样很容易就得到了如下的算法（源码链接）：

func BinarySqrt[T SqrtValue](value T) T {
	if value < 0 {
		panic(fmt.Sprintf("BinarySqrt Negative Value %f", value))
	}
	if value == 0 {
		return 0
	}
	min := T(0.0)
	max := value
	const SquareRootPrecise = 10e-6
	for (max - min) > SquareRootPrecise {
		mid := min + (max-min)/2. // 防治俩数值相加爆掉
		delta := mid*mid - value
		if -ResultPrecise <= delta && delta <= ResultPrecise {
			min = mid
			break
		}
		if delta > 0 {
			max = mid
		} else {
			min = mid
		}
	}
	return min
}

二分法的时间复杂度是log(n)，是一种工程优化后的有效解法，但是它不是一种优质甚至于完美的算法。数学上类似问题，有更好的解决方案——牛顿迭代法。

牛顿迭代法是一种数学上用于实数域和复数域上近似求解方程的方法：只要函数二阶可导，那么在待求的零点x周围存在一个区域，只要起始点x0位于这个邻近区域内，那么迭代必定收敛。

依照牛顿迭代法的定义，求取一个数的平方根，其实就是求f(x) = x ^ 2 - C在f(x)=0(x>=0)时的x的解，选择初始点为x=C (0<= x <= C)（符合初始点要求）。并依据牛顿迭代法，一步步迭代最初的值，即可获得最终满足要求的近似解，下面是一个朴素的牛顿迭代法(源码)。

func NewtonSqrt[T SqrtValue](value T) T {
	if value < 0 {
		panic(fmt.Sprintf("Sqrt A Negative Value %f", value))
	}
	if value == 0 {
		return value
	}
	v := value
	for math.Abs(float64(v*v-value)) > ResultPrecise {
		v = (v + value/v) / 2.
	}
	return v
}

通过后文的benchmark可以发现，朴素的牛顿迭代法的收敛速度要远快于二分法，且它的潜能远远没有到达上限。牛顿迭代法有个很重要的特性，好的起始点的能加速收敛。由此，引出了一段传奇代码——快速平方根倒数算法。看下面那段写在的《雷神之锤III》的引擎里的倒排平方根肃算法（当然不一定是传奇的卡马克写的）：

float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;                       // evil floating point bit level hacking（对浮点数的邪恶位元hack）
	i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//      y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
	return y;
}

代码的核心就是牛顿迭代法，但是它的初始值的计算用了一个神奇的0x5f3759df，通过这个初始值，1～2次迭代就可以得到精度尚可的结果值，它的效率是O(1)。至于为啥这个数字能够获得如此精确的初始值，笔者Google了很多资料，最终在chris lomont的paper找到了答案。当然chris lomont甚至给出了更优的解，笔者“拿来”，并写下了如下代码(源码)：

// QuickSqrt32 倒排平方根
func QuickInvSqrt32(x float32) float32 {
	const MAGIC_NUMBER = 0x5f375a86
	xhalf := 0.5 * x
	i := math.Float32bits(x)
	i = MAGIC_NUMBER - (i >> 1)
	x = math.Float32frombits(i)
	x = x * (1.5 - xhalf*x*x)
	x = x * (1.5 - xhalf*x*x)
	return x
}

// QuickSqrt32 基于快速平方根倒数的开方算法
func QuickSqrt32(x float32) float32 {
	return x * QuickInvSqrt32(x)
}

这么多方法，性能到底各自是个啥级别呢？笔者用benchmark 跑了一下，结果如下：

goos: darwin
goarch: amd64
pkg: go-sqrt/sqrt
cpu: Intel(R) Core(TM) i5-6360U CPU @ 2.00GHz
BenchmarkBinaryFloat32-4         3582926               309.8 ns/op
BenchmarkBinaryFloat64-4         3825048               310.0 ns/op
BenchmarkNewtonFloat32-4         8788845               130.9 ns/op
BenchmarkNewtonFloat64-4         8863969               131.4 ns/op
BenchmarkQuickSqrt32-4          229298521                4.923 ns/op
BenchmarkQuickSqrt64-4          332104941                3.614 ns/op
BenchmarkMath-4                 424857439                2.826 ns/op
BenchmarkMathNotHardware-4       5501313               198.7 ns/op