渐进式优化Geohash算法渐进式优化Ｇｅｏｈａｓｈ算法，最终实现９７．４％的性能提升，提供相关优化过程的理论证明和算法

背景

$\qquad$ Geohash算法是一种地理编码算法，它可以将二维经纬度坐标转换成一维编码，广泛应用于地理坐标存储、搜索。Geohash的编码过程分为三步：

转换经纬度为二进制，首先设定范围 $[start, end]$ ，编码目标为 $target$ 。求解 $mid=\frac{start+end}{2}$ ，如果 $target<mid$ ，则二进制该位置0，同时设 $end=mid$ ；反之，则二进制该位置1，同时设 $start=mid$ 。不断重复这个过程，就得到了经纬度的二进制编码；
按照奇数位放经度二进制，偶数位放纬度二进制的原则，由高到低合并经纬度的二进制编码，生成新的编码；
对新的二进制编码使用Base32算法编码，生成字符串；

以天安门的坐标 $(39.54,116.23)$ 举例，首先对维度进行二进制编码:

start	end	mid	code
-90	90	0	1
0	90	45	0
0	45	22.5	1
22.5	45	33.75	1
33.75	45	39.375	1
39.375	45	42.1875	0
39.375	42.1875	40.78125	0
...	...	...	...

纬度的32位二进制编码为10111000001111000001001100011101。重复上面的过程，得到经度的32位二进制编码为11010010101001110000011010011001。接下来对经度，纬度的二进制编码进行合并，得到二进制编码1110010011000111101111111011111000111110101100001010000001001110。最后进行base32编码，得到9jxzrszc182f。
$\qquad$ Geohash编码属于计算密集型操作，不涉及IO操作，不同语言（Golang, Rust, NodeJS...）的Geohash算法的实现也大体相似。本文启发自 Geohash in Golang Assembly ，相关代码参考了原文的实现。本文的主要工作是对相关算法优化过程提供理论推导。注本文只涉及Geohash算法的1、2过程，不涉及对Base32编码的优化。

二分查找

经纬度编码

func encode(start, end, target float64) uint32 {
   var ans uint32
   for i := 0; i < 32; i++ {
      mid := (start + end) / 2.0
      if target <= mid {
         end = mid
         ans <<= 1
      } else {
         start = mid
         ans = (ans << 1) + 1
      }
   }
   return ans
}

$\qquad$ 对经纬度进行二进制编码可以使用二分查找算法，这是最简单暴力的实现，也是最接近Geohash算法定义的实现。二分查找算法的本质是不停对 $[start,end]$ 进行二分查找，根据 $mid=\frac{start+end}{2}$ 和target大小来设置新的 $start$ 和 $end$ ，最终精准定位 $target$ 在 $[start,end]$ 所在区间。

编码合并

func Spread(x uint32) uint64 {
   X := uint64(x)
   X = (X | (X << 16)) & 0x0000ffff0000ffff
   X = (X | (X << 8)) & 0x00ff00ff00ff00ff
   X = (X | (X << 4)) & 0x0f0f0f0f0f0f0f0f
   X = (X | (X << 2)) & 0x3333333333333333
   X = (X | (X << 1)) & 0x5555555555555555
   return X
}

func Interleave(lat, lng uint32) uint64 {
	return Spread(lat) | (Spread(lng) << 1)
}

$\qquad$ 合并经纬度二进制编码参考了Bit Twiddling Hacks (stanford.edu)。

Benchmark

func EncodeInt(lat, lng float64) uint64 {
   encodeLat := encode(-90.0, 90.0, lat)
   encodeLng := encode(-180.0, 180.0, lng)
   return Interleave(encodeLat, encodeLng)
}

$\qquad$ Benchmark显示基于二分查找的算法实现耗时为 $123.3ns/op$ 。

优化

思考：如何使用二进制表示一个十进制数？

$\qquad$ 在进行算法优化之前，首先思考一个问题：如何使用二进制表示一个十进制数？人类社会是十进制的世界，我们所接触的各种价格，如一台iphone的价格，一个背包的价格均是十进制表示的，而二进制则是计算机世界的基石，计算机中的一切操作都是0|1计算而来。那么一个准确的十进制数是否可以使用二进制数完整表示？为了解答这个问题，可做以下尝试。首先定义 $f=f_{int}+f_{decimal}$ ， $f$ 为浮点数， $f_{int}$ 表示 $f$ 的整数部分，而 $f_{decimal}$ 代表小数部分。整数 $f_{int}$ 可以精确表示为

f_{int}=\sum_0^n{a_i*2^{n-i}},a_i\in\{0,1\}

，即可以使用一个0|1序列表示整数 $f_{int}$ 。具体证明可以使用归纳法，即假设整数 $n$ 可以转换成二进制，那么整数 $n+1$ 必然可以转换成二进制。已知0可以表示成二进制0，由此推导所有整数都可以表示成二进制。
$\qquad$ 整数作为离散的数字，它们之间的差值必然是整数，因此整数绝对可以使用二进制精确表示。而小数作为连续的数字，小数的精度可以精确到小数点后无限位，因此小数肯定无法使用二进制精确表示。但参照整数的二进制表示，可以使用二进制无限逼近，即使用下列等式在“损失精度”的前提下表示一个小数。

f_{decimal}\approx\sum_1^n{a_i*2^{-i}},a_i\in\{0,1\}

以数字0.3举例，如何用二进制表示0.3呢？首先最高位置1，看看表达的小数是否超过0.3。

0.1=1*2^{-1}=0.5

已经超过0.3，遂置0。接下来，次高位继续执行上述操作

0.01=0*2^{-1}+1*2^{-2}=0.25

未超过，次高位置1。不断执行该操作，直至

0.01001=0*2^{-1}+1*2^{-2}+0*2^{-3}+0*2^{-4}+1*2^{-5}=0.28125

不断继续下去，你会发现总是无法精确表示0.3，只能无限逼近0.3。

$\qquad$ 现在尝试用代码实现转换小数为32位二进制，即实现 $f_{decimal}\approx\sum_1^{32}{a_i*2^{-i}},a_i\in\{0,1\}$ ：

func EncodeDecimal(decimal float64) uint32 {
   var ans uint32
   mid := 0.5
   for i := 1; i <= 32; i++ {
      if decimal >= mid {
         ans = (ans << 1) + 1
         decimal -= mid
      } else {
         ans = ans << 1
      }
      mid *= 0.5
   }
   return ans
}

上面的算法使用了for循环，实际上小数的二进制转换还可以加速。尝试做如下等式转换：

f_{decimal}\approx\sum_1^{32}{a_i*2^{-i}},a_i\in\{0,1\} \\ 2^{32}*f_{decimal}\approx\sum_1^{32}{a_i*2^{32-i}},a_i\in\{0,1\} \\ 2^{32}*f_{decimal}\approx\sum_0^{31}{a_i*2^{31-i}},a_i\in\{0,1\}

小数的32位二进制表示等价于将小数放大 $2^{32}$ 倍后的整数的二进制表示。代码实现可以优化成：

func EncodeDecimal2(decimal float64) uint32 {
   return uint32(math.Ldexp(decimal, 32))
}

Benchmark

$\qquad$ 下面是两种算法的Benchmark，通过Benchmark可以发现小小的转换换来了 $87.5\%$ 的性能提升。

二进制编码优化

$\qquad$ 上面绕了这么多，是为了对优化算法的提出进行铺垫。看下面的理论推导，这里经过左右平移、等比缩小等线性运算成功地将经纬度（ $target$ ）的二进制编码转换成求解小数 $\frac{target-start}{end-start}$ 二进制表示。

[start, target, end] \rightrightarrows [0, target-start, end-start] \rightrightarrows [0, \frac{target-start}{end-start},1]

$\qquad$ 经度 $lng$ 的二进制编码等于小数 $\frac{lng}{360.0}+0.5$ 的二进制表示，而纬度 $lat$ 的二进制编码等于 $\frac{lat}{180.0}+0.5$ 的二进制编码。因此，套用上面的求解小数二进制表示的算法实现了下列优化算法：

func EncodeInt2(lat, lng float64) uint64 {
   return Interleave(Quantize(lat, lng))
}

func Quantize(lat, lng float64) (lat32 uint32, lng32 uint32) {
   lat32 = uint32(math.Ldexp((lat+90.0)/180.0, 32))
   lng32 = uint32(math.Ldexp((lng+180.0)/360.0, 32))
   return
}

这里使用了Go的math.Ldexp(float, int)方法。

Benchmark

$\qquad$ 算法的单次计算耗时为 $19.00 ns/op$ ，和前面的基于二分查找的算法相比，性能提升 $84.6\%$ 。

双精度浮点数（IEEE 754）

$\qquad$ Geohash算法的优化仍未完结，上述算法还有优化的空间。本节介绍Geohash算法的标准实现，当然这个实现也是很多语言（Golang, Rust, NodeJS...）的标准实现，标准实现利用了双精度浮点数的表示。双精度浮点数由IEEE 754规范定义，它使用64位二进制表示一个浮点数，64位二进制的拆分如下：

1位	11位	52位
S（符号位）	E（阶码位）	M（小数位）
0表示正，1表示负	1~2046	任意

根据IEEE 754规范，一个浮点数按照以下的公式编码：

f = {(-1)^S}*2^{E-1023}*(1+\sum_1^{52}{M_i*2^{-i}})

优化算法

$\qquad$ 先定义 $y=1+x,x\in[0,1]$ ，即 $y$ 代表一个处于 $[1,2]$ 之间的数。使用IEEE 754规范编码 $y$ ，可进行以下公式推导：

y=1+x \\ y={(-1)^0}*2^{1023-1023}*(1+\sum_1^{52}{M_i*2^{-i}}) \\ y=1+\sum_1^{52}{M_i*2^{-i}} \\ x=\sum_1^{52}{M_i*2^{-i}} \\

$M$ 序列正好是 $x$ 的二进制表示。因此经纬度（ $target$ ）的二进制编码可转换成求 $\frac{target}{2*end}+1.5$ 双精度浮点数的 $M$ 序列。Go中可使用 math.Float64bits(x)获取浮点数x的M序列，基于此进行算法优化：

func QuantizeBits(x, base float64) uint64 {
   return math.Float64bits(x/base + 1.5)
}

func EncodeInt3(lat, lng float64) uint64 {
   return Interleave(uint32(QuantizeBits(lat, 180.0)>>20), uint32(QuantizeBits(lat, 360.0)>>20))
}

Benchmark

$\qquad$ Oops! 算法耗时由 $19ns/op$ 下降到 $7.73ns/op$ ，性能提升 $59.3\%$ 。

汇编重构

$\qquad$ 当算法优化到一种程度的时候，可以考虑使用汇编重构算法，来进一步提升算法性能。首先进行宏定义：

// +build amd64,go1.6

#include "textflag.h"

#define LATF   X0
#define LATI   R8
#define LNGF   X1
#define LNGI   R9
#define TEMP   R10
#define GHSH   R11
#define MASK   BX

重构

二进制编码重构

// func QuantizeAsm(lat, lng float64) (uint32, uint32)
TEXT ·QuantizeAsm(SB), NOSPLIT, $0

    MOVSD lat+0(FP), LATF
    MOVSD lng+8(FP), LNGF


    MULSD $(0.005555555555555556), LATF
    ADDSD $(1.5), LATF                 // lat/180.0+1.5

    MULSD $(0.002777777777777778), LNGF
    ADDSD $(1.5), LNGF                // lng/360.0+1.5

    MOVQ LNGF, LNGI
    SHRQ $20, LNGI                   // 右移20位

    MOVQ  LATF, LATI
    SHRQ  $20, LATI                 // 右移20位

    MOVQ LATI, ret+16(FP)
    MOVQ LNGI, ret+20(FP)
    RET

编码合并重构

// func InterleaveAsm(x, y uint32) uint64
TEXT ·InterleaveAsm(SB), NOSPLIT, $0


    MOVSD x+0(FP), LATF
    MOVSD x+4(FP), LNGF


    MOVQ $0x5555555555555555, MASK  // 010101010101010101010101010101010101(64位)
    MOVQ LATF, LATI
    MOVQ LNGF, LNGI

    PDEPQ MASK, LATI, GHSH         // lat二进制编码按照MASK映射
    PDEPQ MASK, LNGI, TEMP         // lng二进制编码按照MASK映射

    SHLQ $1, TEMP                  // lng映射后编码左移一位
    XORQ TEMP, GHSH                // 或操作进行编码合并

    MOVQ GHSH, ret+8(FP)
    RET

合并


// func EncodeIntAsm(lat, lng float64) uint64
TEXT ·EncodeIntAsm(SB), NOSPLIT, $0
   CMPB ·useAsm(SB), $1
   JNE  fallback

   MOVSD lat+0(FP), LATF
   MOVSD lng+8(FP), LNGF

   MOVQ $0x5555555555555555, MASK

   MULSD $(0.005555555555555556), LATF
   ADDSD $(1.5), LATF

   MULSD $(0.002777777777777778), LNGF
   ADDSD $(1.5), LNGF

   MOVQ LNGF, LNGI
   SHRQ $20, LNGI

   MOVQ  LATF, LATI
   SHRQ  $20, LATI
   PDEPQ MASK, LATI, GHSH

   PDEPQ MASK, LNGI, TEMP

   SHLQ $1, TEMP
   XORQ TEMP, GHSH

   MOVQ GHSH, ret+16(FP)
   RET

fallback:
   JMP ·encodeInt(SB)

Benchmark

$\qquad$ 通过Benchmark跑分发现，算法经过汇编重构后，耗时由 $7.73 ns/op$ 下降到了 $2.82 ns/op$ ，性能提升了 $63.5\%$ 。

总结

$\qquad$ Geohash算法经过一系列的优化，单位操作耗时由 $123.7ns/op$ 下降到 $2.794ns/op$ ，性能提升 $97.74\%$ ，优化效果十分明显。Geohash算法的优化思路说明，将一个问题转换成另一个问题可能会收获意想不到的效果，很多“伟大”的优化本质上是另一维度的“投机取巧”，算法优化的终极有可能是汇编。

一种Rust实现

$\qquad$ 这里提供一种Rust版本的Geohash实现，当然你也可以试着用其他语言参照下面的算法实现。

// 扩散
pub fn spread(x: u32) -> u64 {
    let mut ans = x as u64;
    ans = (ans | (ans << 16)) & 0x0000ffff0000ffff;
    ans = (ans | (ans << 8)) & 0x00ff00ff00ff00ff;
    ans = (ans | (ans << 4)) & 0x0f0f0f0f0f0f0f0f;
    ans = (ans | (ans << 2)) & 0x3333333333333333;
    ans = (ans | (ans << 1)) & 0x5555555555555555;
    ans
}

// 合并
pub fn interleave(lat: u32, lng: u32) -> u64 {
    spread(lat) | (spread(lng) << 1)
}

// 归一化
pub fn quantize(lat: f64, base: f64) -> u32 {
    ((lat/base + 1.5).to_bits() >> 20) as u32
}



pub fn encode(lat: f64, lng: f64) -> u64 {
    interleave(quantize(lat, 180.0), quantize(lng, 360.0))
}

Benchmark

$\qquad$ Benchmark显示Rust版本的Geohash算法耗时为 $7 ns/op$ 。

渐进式优化Geohash算法

背景

二分查找

经纬度编码

编码合并

Benchmark

优化

思考：如何使用二进制表示一个十进制数？

Benchmark

二进制编码优化

Benchmark

双精度浮点数（IEEE 754）

优化算法

Benchmark

汇编重构

重构

二进制编码重构

编码合并重构

合并

Benchmark

总结

一种Rust实现

Benchmark

参考文献