Speedup

dbms/src/Functions/greatCircleDistance.cpp

@@ -26,6 +26,91 @@ namespace ErrorCodes
 
 static inline Float64 degToRad(Float64 angle) { return angle * DEGREES_IN_RADIANS; }
 
+static const double PI = 3.14159265358979323846;
+static const double TO_RAD = PI / 180.0;
+static const double TO_RAD2 = PI / 360.0;
+static const double TO_DEG = 180.0 / PI;
+static const float TO_RADF = static_cast<float>(PI / 180.0);
+static const float TO_RADF2 = static_cast<float>(PI / 360.0);
+static const float TO_DEGF = static_cast<float>(180.0 / PI);
+
+const int GEODIST_TABLE_COS = 1024; // maxerr 0.00063%
+const int GEODIST_TABLE_ASIN = 512;
+const int GEODIST_TABLE_K = 1024;
+
+static float g_GeoCos[GEODIST_TABLE_COS + 1];        ///< cos(x) table
+static float g_GeoAsin[GEODIST_TABLE_ASIN + 1];    ///< asin(sqrt(x)) table
+static float g_GeoFlatK[GEODIST_TABLE_K + 1][2];    ///< GeodistAdaptive() flat ellipsoid method k1,k2 coeffs table
+
+inline double sqr ( double v ) { return v*v;}
+inline float fsqr ( float v ) { return v*v; }
+
+
+void GeodistInit()
+{
+    for (int i = 0; i <= GEODIST_TABLE_COS; i++)
+        g_GeoCos[i] = (float) cos(2 * PI * i / GEODIST_TABLE_COS); // [0, 2pi] -> [0, COSTABLE]
+
+    for (int i = 0; i <= GEODIST_TABLE_ASIN; i++)
+        g_GeoAsin[i] = (float) asin(sqrt(double(i) / GEODIST_TABLE_ASIN)); // [0, 1] -> [0, ASINTABLE]
+
+    for (int i = 0; i <= GEODIST_TABLE_K; i++)
+    {
+        double x = PI * i / GEODIST_TABLE_K - PI * 0.5; // [-pi/2, pi/2] -> [0, KTABLE]
+        g_GeoFlatK[i][0] = (float) sqr(111132.09 - 566.05 * cos(2 * x) + 1.20 * cos(4 * x));
+        g_GeoFlatK[i][1] = (float) sqr(111415.13 * cos(x) - 94.55 * cos(3 * x) + 0.12 * cos(5 * x));
+    }
+}
+
+static inline float GeodistDegDiff(float f)
+{
+    f = (float) fabs(f);
+    while (f > 360)
+        f -= 360;
+    if (f > 180)
+        f = 360 - f;
+    return f;
+}
+
+static inline float GeodistFastCos(float x)
+{
+    float y = (float) (fabs(x) * GEODIST_TABLE_COS / PI / 2);
+    int i = int(y);
+    y -= i;
+    i &= (GEODIST_TABLE_COS - 1);
+    return g_GeoCos[i] + (g_GeoCos[i + 1] - g_GeoCos[i]) * y;
+}
+
+static inline float GeodistFastSin(float x)
+{
+    float y = float(fabs(x) * GEODIST_TABLE_COS / PI / 2);
+    int i = int(y);
+    y -= i;
+    i = (i - GEODIST_TABLE_COS / 4) & (GEODIST_TABLE_COS - 1); // cos(x-pi/2)=sin(x), costable/4=pi/2
+    return g_GeoCos[i] + (g_GeoCos[i + 1] - g_GeoCos[i]) * y;
+}
+
+
+/// fast implementation of asin(sqrt(x))
+/// max error in floats 0.00369%, in doubles 0.00072%
+static inline float GeodistFastAsinSqrt(float x)
+{
+    if (x < 0.122)
+    {
+        // distance under 4546km, Taylor error under 0.00072%
+        float y = (float) sqrt(x);
+        return y + x * y * 0.166666666666666f + x * x * y * 0.075f + x * x * x * y * 0.044642857142857f;
+    }
+    if (x < 0.948)
+    {
+        // distance under 17083km, 512-entry LUT error under 0.00072%
+        x *= GEODIST_TABLE_ASIN;
+        int i = int(x);
+        return g_GeoAsin[i] + (g_GeoAsin[i + 1] - g_GeoAsin[i]) * (x - i);
+    }
+    return (float) asin(sqrt(x)); // distance over 17083km, just compute honestly
+}
+
 /**
  *  The function calculates distance in meters between two points on Earth specified by longitude and latitude in degrees.
  *  The function uses great circle distance formula https://en.wikipedia.org/wiki/Great-circle_distance.
@@ -38,7 +123,7 @@ class FunctionGreatCircleDistance : public IFunction
 public:
 
     static constexpr auto name = "greatCircleDistance";
-    static FunctionPtr create(const Context &) { return std::make_shared<FunctionGreatCircleDistance>(); }
+    static FunctionPtr create(const Context &) { GeodistInit(); return std::make_shared<FunctionGreatCircleDistance>(); }
 
 private:
 
@@ -103,16 +188,38 @@ private:
             lat1Deg < -90 || lat1Deg > 90 ||
             lat2Deg < -90 || lat2Deg > 90)
         {
-            throw Exception("Arguments values out of bounds for function " + getName(), ErrorCodes::ARGUMENT_OUT_OF_BOUND);
+            throw Exception("Arguments values out of bounds for function " + getName(),
+                            ErrorCodes::ARGUMENT_OUT_OF_BOUND);
         }
 
-        Float64 lon1Rad = degToRad(lon1Deg);
+        float dlat = GeodistDegDiff(lat1Deg - lat2Deg);
-        Float64 lat1Rad = degToRad(lat1Deg);
+        float dlon = GeodistDegDiff(lon1Deg - lon2Deg);
-        Float64 lon2Rad = degToRad(lon2Deg);
+
-        Float64 lat2Rad = degToRad(lat2Deg);
+        if (dlon < 13)
-        Float64 u = sin((lat2Rad - lat1Rad) / 2);
+        {
-        Float64 v = sin((lon2Rad - lon1Rad) / 2);
+            // points are close enough; use flat ellipsoid model
-        return 2.0 * EARTH_RADIUS_IN_METERS * asin(sqrt(u * u + cos(lat1Rad) * cos(lat2Rad) * v * v));
+            // interpolate sqr(k1), sqr(k2) coefficients using latitudes midpoint
+            float m = (lat1Deg + lat2Deg + 180) * GEODIST_TABLE_K / 360; // [-90, 90] degrees -> [0, KTABLE] indexes
+            int i = int(m);
+            i &= (GEODIST_TABLE_K - 1);
+            float kk1 = g_GeoFlatK[i][0] + (g_GeoFlatK[i + 1][0] - g_GeoFlatK[i][0]) * (m - i);
+            float kk2 = g_GeoFlatK[i][1] + (g_GeoFlatK[i + 1][1] - g_GeoFlatK[i][1]) * (m - i);
+            return (float) sqrt(kk1 * dlat * dlat + kk2 * dlon * dlon);
+        }
+        // points too far away; use haversine
+        static const float D = 2 * 6371000;
+        float a = fsqr(GeodistFastSin(dlat * TO_RADF2)) +
+                  GeodistFastCos(lat1Deg * TO_RADF) * GeodistFastCos(lat2Deg * TO_RADF) *
+                  fsqr(GeodistFastSin(dlon * TO_RADF2));
+        return (float) (D * GeodistFastAsinSqrt(a));
+
+//        Float64 lon1Rad = degToRad(lon1Deg);
+//        Float64 lat1Rad = degToRad(lat1Deg);
+//        Float64 lon2Rad = degToRad(lon2Deg);
+//        Float64 lat2Rad = degToRad(lat2Deg);
+//        Float64 u = sin((lat2Rad - lat1Rad) / 2);
+//        Float64 v = sin((lon2Rad - lon1Rad) / 2);
+//        return 2.0 * EARTH_RADIUS_IN_METERS * asin(sqrt(u * u + cos(lat1Rad) * cos(lat2Rad) * v * v));
     }