Added uint64_t geohash implementation with float black magic

Author: Lezh1k

C/include/Geohash.h

@@ -4,17 +4,12 @@
 #include <stdint.h>
 
 static const int GEOHASH_MAX_PRECISION = 12;
-void GeohashEncode(double srcLat,
-                   double srcLon,
-                   char* geohash,
-                   uint8_t precision);
+
+uint64_t GeohashEncodeU64(double lat, double lon);
+int GeohashComparePointsU64(double lon1, double lat1, double lon2, double lat2, int precision);
 
 void GeohashDecode(const char* str,
                    double *pLon,
                    double *pLat);
 
-int GeohashComparePoints(double lon1, double lat1,
-                         double lon2, double lat2,
-                         int precision);
-
 #endif // GEOHASH_H

C/src/Coordinates.cpp

@@ -44,23 +44,23 @@ CoordFilterByGeoHash(std::vector<geopoint_t> &lstSrc,
     int32_t index;
     int32_t count;
   };
-  static char buff[GEOHASH_MAX_PRECISION+1] = {0};
 
   std::vector<geopoint_t> lstRes;
-  std::map<std::string, AuxItem> dctHashCount;
-  typedef std::map<std::string, AuxItem>::iterator dctIter;
+  std::map<uint64_t, AuxItem> dctHashCount;
+  typedef std::map<uint64_t, AuxItem>::iterator dctIter;
 
   int idx = 0;
   for (auto ci = lstSrc.begin(); ci != lstSrc.end(); ++ci) {
-    GeohashEncode(ci->Latitude, ci->Longitude, buff, precision);
-    std::string geohash(buff, precision);
-    dctIter it = dctHashCount.find(geohash);
+    uint64_t gh = GeohashEncodeU64(ci->Latitude, ci->Longitude);
+    gh >>= precision*5 + 4;
+
+    dctIter it = dctHashCount.find(gh);
     if (it == dctHashCount.end()) {
       AuxItem ni;
       ni.count = 0;
       ni.lat = ni.lon = 0.0;
       ni.index = NOT_VALID_POINT_INDEX;
-      auto ir = dctHashCount.insert(std::pair<std::string, AuxItem>(geohash, ni));
+      auto ir = dctHashCount.insert(std::pair<uint64_t, AuxItem>(gh, ni));
       if (!ir.second)
         continue;
       it = ir.first;

C/src/Geohash.cpp

@@ -4,113 +4,50 @@
 #include <string.h>
 #include <stdio.h>
 
+#include <iostream>
 #include <immintrin.h>
 #include "Geohash.h"
 
-//where is 'A'???? oh my god, we lost 'I', 'L', 'O' too. :(
-#define BASE32_COUNT 32
-static const char base32Table[BASE32_COUNT] = {'0', '1', '2', '3', '4', '5', '6', '7',
-                                               '8', '9', 'b', 'c', 'd', 'e', 'f', 'g',
-                                               'h', 'j', 'k', 'm', 'n', 'p', 'q', 'r',
-                                               's', 't', 'u', 'v', 'w', 'x', 'y', 'z'};
-static int searchInBase32Table(char c) {
-  int f, l, m;
-  f = 0;
-  l = BASE32_COUNT;
-  while (f < l) {
-    m = (f + l) >> 1;
-    if (c <= base32Table[m])
-      l = m;
-    else
-      f = m+1;
-  }
-  return l;
-}
-//////////////////////////////////////////////////////////////////////////
-
-typedef struct interval {
-  double min, max;
-} interval_t;
-///////////////////////////////////////////////////////
-
-void GeohashEncode(double srcLat,
-                   double srcLon,
-                   char* geohash,
-                   uint8_t precision) {
-  interval_t lat = {-90.0, 90.0};
-  interval_t lon = {-180.0, 180.0};
-  interval_t *ci;
-  bool isEven = true;
-  double mid, cd;
-  int32_t idx = 0; // index into base32 map
-  int8_t bit = 0; // each char holds 5 bits
+static uint64_t interleave(uint64_t x, uint64_t y) {
+  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;
 
-  while (precision) {
-    if (isEven) {
-      ci = &lon;
-      cd = srcLon;
-    } else {
-      ci = &lat;
-      cd = srcLat;
-    }
+  y = (y | (y << 16)) & 0x0000ffff0000ffff;
+  y = (y | (y << 8)) & 0x00ff00ff00ff00ff;
+  y = (y | (y << 4)) & 0x0f0f0f0f0f0f0f0f;
+  y = (y | (y << 2)) & 0x3333333333333333;
+  y = (y | (y << 1)) & 0x5555555555555555;
 
-    mid = (ci->min + ci->max) / 2.0;
-    idx <<= 1; //idx *= 2
-    if (cd >= mid) {
-      ci->min = mid;
-      idx += 1;
-    } else {
-      ci->max = mid;
-    }
+  return x | (y << 1);
 
-    isEven = !isEven;
-
-    if (++bit == 5) {
-      *geohash++ = base32Table[idx];
-      idx = bit = 0;
-      --precision;
-    }
-  }
-  *geohash++ = 0;
+  //use pdep instructions
+//  return _pdep_u64(x, 0x5555555555555555) | _pdep_u64(y, 0xaaaaaaaaaaaaaaaa);
 }
-//////////////////////////////////////////////////////////////////////////
-
-void GeohashDecode(const char* str,
-                   double *pLon,
-                   double *pLat) {
-  interval_t latInterval = {-90.0, 90.0};
-  interval_t lonInterval = {-180.0, 180.0};
-  interval_t *ci;
-  bool isEven = true;
-  int idx ;
-  double mid;
+///////////////////////////////////////////////////////
 
-  for (; *str; ++str) {
-    idx = searchInBase32Table(*str);
-    for (int i = 0; i < 5; ++i) {
-      ci = isEven ? &lonInterval : &latInterval;
-      mid = (ci->max - ci->min) / 2.0;
-      if ((idx << i) & 0x10) {
-        ci->min += mid;
-      } else {
-        ci->max -= mid;
-      }
-      isEven = !isEven;
-    }
-  } //for c in str
-  *pLat = latInterval.max - ((latInterval.max - latInterval.min) / 2.0);
-  *pLon = lonInterval.max - ((lonInterval.max - lonInterval.min) / 2.0);
+uint64_t GeohashEncodeU64(double lat, double lon) {
+  lat = lat/180.0 + 1.5;
+  lon = lon/360.0 + 1.5;
+  uint64_t ilat = *((uint64_t*)&lat);
+  uint64_t ilon = *((uint64_t*)&lon);
+  ilat >>= 20;
+  ilon >>= 20;
+  ilat &= 0x00000000ffffffff;
+  ilon &= 0x00000000ffffffff;
+
+  return interleave(ilat, ilon);
 }
-//////////////////////////////////////////////////////////////////////////
 
-int GeohashComparePoints(double lon1, double lat1,
-                         double lon2, double lat2,
-                         int precision) {
+int GeohashComparePointsU64(double lon1, double lat1,
+                            double lon2, double lat2,
+                            int precision) {
   assert(precision >= 1 && precision <= GEOHASH_MAX_PRECISION);
-  static char geohash1[GEOHASH_MAX_PRECISION+1] = {0};
-  static char geohash2[GEOHASH_MAX_PRECISION+1] = {0};
-
-  GeohashEncode(lat1, lon1, geohash1, precision);
-  GeohashEncode(lat2, lon2, geohash2, precision);
-  return memcmp(geohash1, geohash2, precision);
+  uint64_t gh1 = GeohashEncodeU64(lat1, lon1);
+  uint64_t gh2 = GeohashEncodeU64(lat2, lon2);
+  gh1 >>= precision*5 + 4;
+  gh2 >>= precision*5 + 4;
+  return gh1 - gh2;
 }

C/src/mainwindow.cpp

@@ -101,33 +101,34 @@ static QString jsCoordsString(const std::vector<geopoint_t>& lst,
 static double filterDistanceRealTime(const std::vector<geopoint_t> &lst,
                                      int prec,
                                      int minPointCount) {
+  static const double COORD_NOT_INITIALIZED = 361.0;
+
   if (lst.empty() || lst.size() == 1)
     return 0.0;  
 
   double distance = 0.0;
-  char buff1[GEOHASH_MAX_PRECISION+1] = {0};
-  char buff2[GEOHASH_MAX_PRECISION+1] = {0};
-  char *pbuff1 = buff1;
-  char *pbuff2 = buff2;
-  char **ppReadGeohash = &pbuff1;
-  char **ppCompGeohash = &pbuff2;
-
   int count;
   geopoint_t tmpGeo, laGeo;
   auto pi = lst.begin();
+  uint64_t gh0, gh;
+  uint64_t *tgh0, *tgh;
+  tgh0 = &gh0;
+  tgh = &gh;
 
-  static const double COORD_NOT_INITIALIZED = 361.0;
   laGeo.Latitude = laGeo.Longitude = COORD_NOT_INITIALIZED;
-  GeohashEncode(pi->Latitude, pi->Longitude, *ppCompGeohash, prec);
+
+  *tgh0 = GeohashEncodeU64(pi->Latitude, pi->Longitude);
   tmpGeo.Latitude = pi->Latitude;
   tmpGeo.Longitude = pi->Longitude;
   count = 1;
 
+  *tgh >>= prec*5 + 4;
   for (++pi; pi != lst.end(); ++pi) {
-    GeohashEncode(pi->Latitude, pi->Longitude, *ppReadGeohash, prec);
+     *tgh = GeohashEncodeU64(pi->Latitude, pi->Longitude);
+     *tgh >>= prec*5 + 4;
 
     //if (ppCompGeohash != ppReadGeohash)
-    if (memcmp(*ppCompGeohash, *ppReadGeohash, prec)) {
+    if (*tgh - *tgh0) {
       if (count >= minPointCount) {
         tmpGeo.Latitude /= count;
         tmpGeo.Longitude /= count;
@@ -144,7 +145,7 @@ static double filterDistanceRealTime(const std::vector<geopoint_t> &lst,
       count = 1;
       tmpGeo.Latitude = pi->Latitude;
       tmpGeo.Longitude = pi->Longitude;
-      std::swap(*ppCompGeohash, *ppReadGeohash);
+      std::swap(*tgh0, *tgh);
       continue;
     }