From dec4701ec727de78728c839ed91fd54cf30337f3 Mon Sep 17 00:00:00 2001 From: Steven Braeger Date: Mon, 15 Sep 2014 14:10:26 -0400 Subject: [PATCH 1/4] Update geohash.h The API has no reason to pass the range arguments by reference. Doing so slows down the function. The arguments are not return values, and the first thing the code does is dereference them and copy them into stack variables. Passing the input ranges by value skips this step and has fewer errors. Furthermore, the arguments are small enough (2 primatives) that putting them onto the stack by value can make the function call inline even faster. --- geohash.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/geohash.h b/geohash.h index 973db8b..d76842f 100755 --- a/geohash.h +++ b/geohash.h @@ -84,11 +84,11 @@ extern "C" * 0:success * -1:failed */ - int geohash_encode(const GeoHashRange* lat_range, const GeoHashRange* lon_range, double latitude, double longitude, uint8_t step, GeoHashBits* hash); - int geohash_decode(const GeoHashRange* lat_range, const GeoHashRange* lon_range, const GeoHashBits* hash, GeoHashArea* area); - int geohash_get_neighbors(const GeoHashBits* hash, GeoHashNeighbors* neighbors); + int geohash_encode(GeoHashRange lat_range,GeoHashRange lon_range, double latitude, double longitude, uint8_t step, GeoHashBits* hash); + int geohash_decode(GeoHashRange lat_range,GeoHashRange lon_range, GeoHashBits hash, GeoHashArea* area); + int geohash_get_neighbors(GeoHashBits hash, GeoHashNeighbors* neighbors); - int geohash_get_neighbor(const GeoHashBits* hash, GeoDirection direction, GeoHashBits* neighbor); + int geohash_get_neighbor(GeoHashBits hash, GeoDirection direction, GeoHashBits* neighbor); GeoHashBits geohash_next_leftbottom(GeoHashBits bits); GeoHashBits geohash_next_rightbottom(GeoHashBits bits); From d17198588076086f91a6f86535b1fbffbb767c3a Mon Sep 17 00:00:00 2001 From: Steven Braeger Date: Mon, 15 Sep 2014 14:13:41 -0400 Subject: [PATCH 2/4] Update geohash.c Same change as before: change input arguments to stack based arguments for efficiency. This also allows us to skip the case where the inputs == NULL in the error check --- geohash.c | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/geohash.c b/geohash.c index 41614b1..6ef0b44 100755 --- a/geohash.c +++ b/geohash.c @@ -44,16 +44,13 @@ * ----------------- */ -int geohash_encode(const GeoHashRange* lat_r, const GeoHashRange* lon_r, double latitude, double longitude, +int geohash_encode(GeoHashRange lat_range, GeoHashRange lon_range, double latitude, double longitude, uint8_t step, GeoHashBits* hash) { - if (NULL == hash || step > 32 || step == 0 || NULL == lat_r || NULL == lon_r) + if (NULL == hash || step > 32 || step == 0) { return -1; } - GeoHashRange lat_range, lon_range; - lat_range = *lat_r; - lon_range = *lon_r; hash->bits = 0; hash->step = step; uint8_t i = 0; @@ -98,19 +95,19 @@ static inline uint8_t get_bit(uint64_t bits, uint8_t pos) return (bits >> pos) & 0x01; } -int geohash_decode(const GeoHashRange* lat_range, const GeoHashRange* lon_range, const GeoHashBits* hash, +int geohash_decode(GeoHashRange lat_range, GeoHashRange lon_range, const GeoHashBits* hash, GeoHashArea* area) { - if (NULL == hash || NULL == area || NULL == lat_range || NULL == lon_range) + if (NULL == hash || NULL == area) { return -1; } area->hash = *hash; uint8_t i = 0; - area->latitude.min = lat_range->min; - area->latitude.max = lat_range->max; - area->longitude.min = lon_range->min; - area->longitude.max = lon_range->max; + area->latitude.min = lat_range.min; + area->latitude.max = lat_range.max; + area->longitude.min = lon_range.min; + area->longitude.max = lon_range.max; for (; i < hash->step; i++) { uint8_t lat_bit, lon_bit; @@ -180,7 +177,7 @@ static int geohash_move_y(GeoHashBits* hash, int8_t d) return 0; } -int geohash_get_neighbors(const GeoHashBits* hash, GeoHashNeighbors* neighbors) +int geohash_get_neighbors(GeoHashBits hash, GeoHashNeighbors* neighbors) { geohash_get_neighbor(hash, GEOHASH_NORTH, &neighbors->north); geohash_get_neighbor(hash, GEOHASH_EAST, &neighbors->east); @@ -193,13 +190,13 @@ int geohash_get_neighbors(const GeoHashBits* hash, GeoHashNeighbors* neighbors) return 0; } -int geohash_get_neighbor(const GeoHashBits* hash, GeoDirection direction, GeoHashBits* neighbor) +int geohash_get_neighbor(GeoHashBits hash, GeoDirection direction, GeoHashBits* neighbor) { if (NULL == neighbor) { return -1; } - *neighbor = *hash; + *neighbor = hash; switch (direction) { case GEOHASH_NORTH: From ec1c7e0e0162d956bf8dd969c3b35388fbe0e846 Mon Sep 17 00:00:00 2001 From: Steven Braeger Date: Mon, 15 Sep 2014 14:57:54 -0400 Subject: [PATCH 3/4] Update geohash.c --- geohash.c | 157 ++++++++++++++++++++++++++++++++++-------------------- 1 file changed, 98 insertions(+), 59 deletions(-) diff --git a/geohash.c b/geohash.c index 6ef0b44..225445e 100755 --- a/geohash.c +++ b/geohash.c @@ -44,6 +44,63 @@ * ----------------- */ +static inline uint64_t interleave64(uint32_t xlo,uint32_t ylo) +{ + static const uint64_t B[] = {0x5555555555555555, 0x3333333333333333, 0x0F0F0F0F0F0F0F0F, 0x00FF00FF00FF00FF, 0x0000FFFF0000FFFF}; + static const unsigned int S[] = {1, 2, 4, 8, 16}; + + uint64_t x=xlo; // Interleave lower bits of x and y, so the bits of x + uint64_t y=ylo; // are in the even positions and bits from y in the odd; //https://graphics.stanford.edu/~seander/bithacks.html#InterleaveBMN + + // x and y must initially be less than 2**32. + + x = (x | (x << S[4])) & B[4]; + y = (y | (y << S[4])) & B[4]; + + x = (x | (x << S[3])) & B[3]; + y = (y | (y << S[3])) & B[3]; + + x = (x | (x << S[2])) & B[2]; + y = (y | (y << S[2])) & B[2]; + + x = (x | (x << S[1])) & B[1]; + y = (y | (y << S[1])) & B[1]; + + x = (x | (x << S[0])) & B[0]; + y = (y | (y << S[0])) & B[0]; + + return x | (y << 1); +} + +static inline uint64_t deinterleave64(uint64_t interleaved) +{ + static const uint64_t B[] = {0x5555555555555555, 0x3333333333333333, 0x0F0F0F0F0F0F0F0F, 0x00FF00FF00FF00FF, 0x0000FFFF0000FFFF,0x00000000FFFFFFFF}; + static const unsigned int S[] = {0, 1, 2, 4, 8, 16}; + + uint64_t x=interleaved; ///reverse the interleave process (http://stackoverflow.com/questions/4909263/how-to-efficiently-de-interleave-bits-inverse-morton) + uint64_t y=interleaved; + + x = (x | (x >> S[0])) & B[0]; + y = (y | (y >> S[0])) & B[0]; + + x = (x | (x >> S[1])) & B[1]; + y = (y | (y >> S[1])) & B[1]; + + x = (x | (x >> S[2])) & B[2]; + y = (y | (y >> S[2])) & B[2]; + + x = (x | (x >> S[3])) & B[3]; + y = (y | (y >> S[3])) & B[3]; + + x = (x | (x >> S[4])) & B[4]; + y = (y | (y >> S[4])) & B[4]; + + x = (x | (x >> S[5])) & B[5]; + y = (y | (y >> S[5])) & B[5]; + + return x | (y << 32); +} + int geohash_encode(GeoHashRange lat_range, GeoHashRange lon_range, double latitude, double longitude, uint8_t step, GeoHashBits* hash) { @@ -58,35 +115,24 @@ int geohash_encode(GeoHashRange lat_range, GeoHashRange lon_range, double latitu { return -1; } - - for (; i < step; i++) - { - uint8_t lat_bit, lon_bit; - if (lat_range.max - latitude >= latitude - lat_range.min) - { - lat_bit = 0; - lat_range.max = (lat_range.max + lat_range.min) / 2; - } - else - { - lat_bit = 1; - lat_range.min = (lat_range.max + lat_range.min) / 2; - } - if (lon_range.max - longitude >= longitude - lon_range.min) - { - lon_bit = 0; - lon_range.max = (lon_range.max + lon_range.min) / 2; - } - else - { - lon_bit = 1; - lon_range.min = (lon_range.max + lon_range.min) / 2; - } - hash->bits <<= 1; - hash->bits += lon_bit; - hash->bits <<= 1; - hash->bits += lat_bit; - } + + + //the algorithm computes the morton code for the geohash location within the range. + //this can be done MUCH more efficiently using the following code + + //compute the coordinate in the range 0-1 + double lat_offset=(latitude-lat_range.min)/(lat_range.max-lat_range.min); + double lon_offset=(longitude-lon_range.min)/(lon_range.max-lon_range.min); + + //convert it to fixed point based on the step size + lat_offset *= (1 << step); + lon_offset *= (1 << step); + + uint32_t ilato=(uint32_t)lat_offset; + uint32_t ilono=(uint32_t)lon_offset; + + //interleave the bits to create the morton code. No branching and no bounding + hash->bits = interleave64(ilato,ilono); return 0; } @@ -95,41 +141,34 @@ static inline uint8_t get_bit(uint64_t bits, uint8_t pos) return (bits >> pos) & 0x01; } -int geohash_decode(GeoHashRange lat_range, GeoHashRange lon_range, const GeoHashBits* hash, +int geohash_decode(GeoHashRange lat_range, GeoHashRange lon_range, GeoHashBits hash, GeoHashArea* area) { - if (NULL == hash || NULL == area) + if (NULL == area) { return -1; } - area->hash = *hash; - uint8_t i = 0; - area->latitude.min = lat_range.min; - area->latitude.max = lat_range.max; - area->longitude.min = lon_range.min; - area->longitude.max = lon_range.max; - for (; i < hash->step; i++) - { - uint8_t lat_bit, lon_bit; - lon_bit = get_bit(hash->bits, (hash->step - i) * 2 - 1); - lat_bit = get_bit(hash->bits, (hash->step - i) * 2 - 2); - if (lat_bit == 0) - { - area->latitude.max = (area->latitude.max + area->latitude.min) / 2; - } - else - { - area->latitude.min = (area->latitude.max + area->latitude.min) / 2; - } - if (lon_bit == 0) - { - area->longitude.max = (area->longitude.max + area->longitude.min) / 2; - } - else - { - area->longitude.min = (area->longitude.max + area->longitude.min) / 2; - } - } + area->hash = hash; + uint8_t step=hash.step; + uint64_t xyhilo=deinterleave64(hash.bits); //decode morton code + + double lat_scale=lat_range.max-lat_range.min; + double lon_scale=lon_range.max-lon_range.min; + + ilato=xyhilo; //get back the original integer coordinates + ilono=xyhilo >> 32; + + //double lat_offset=ilato; + //double lon_offset=ilono; + //lat_offset /= (1<latitude.min = lat_range.min+ldexp(ilato,-step)*lat_scale; + area->latitude.max = lat_range.min+ldexp(ilato+1,-step)*lat_scale; + area->longitude.min = lon_range.min+ldexp(ilono,-step)*lon_scale; + area->longitude.max = lon_range.min+ldexp(ilono+1,-step)*lon_scale; + return 0; } From 768f561754de4d0da6b0a3a21f96ecb96325bd96 Mon Sep 17 00:00:00 2001 From: Steven Braeger Date: Mon, 15 Sep 2014 14:59:20 -0400 Subject: [PATCH 4/4] Update geohash.c I re-wrote the encode and decode functions to provide a DRAMATICALLY faster interface that does no branching and uses only a few 64-bit integer operations with unrolled loops to do the interleaving/deinterleaving using binary magic numbers for efficient morton codes. The code should go much much much faster now. --- geohash.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geohash.c b/geohash.c index 225445e..f336da8 100755 --- a/geohash.c +++ b/geohash.c @@ -43,7 +43,7 @@ * | 0,0 | 1,0 | * ----------------- */ - +//updated files. static inline uint64_t interleave64(uint32_t xlo,uint32_t ylo) { static const uint64_t B[] = {0x5555555555555555, 0x3333333333333333, 0x0F0F0F0F0F0F0F0F, 0x00FF00FF00FF00FF, 0x0000FFFF0000FFFF};