From e48ac075c053c39d3f7ed4d86deb6554b7d97e1a Mon Sep 17 00:00:00 2001
From: filipe oliveira <filipecosta.90@gmail.com>
Date: Mon, 5 Dec 2022 13:45:04 +0000
Subject: [PATCH] GEOSEARCH BYBOX: Simplified haversine distance formula when
 longitude diff is 0 (#11579)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

This is take 2 of `GEOSEARCH BYBOX` optimizations based on haversine
distance formula when longitude diff is 0.
The first one was in #11535 .

- Given longitude diff is 0 the asin(sqrt(a)) on the haversine is asin(sin(abs(u))).
- arcsin(sin(x)) equal to x when x ∈[−𝜋/2,𝜋/2].
- Given latitude is between [−𝜋/2,𝜋/2] we can simplifiy arcsin(sin(x)) to x.

On the sample dataset with 60M datapoints, we've measured 55% increase
in the achievable ops/sec.
---
 src/geohash_helper.c | 28 ++++++++++++++++++----------
 1 file changed, 18 insertions(+), 10 deletions(-)

diff --git a/src/geohash_helper.c b/src/geohash_helper.c
index a089c97c2..a3816fbe3 100644
--- a/src/geohash_helper.c
+++ b/src/geohash_helper.c
@@ -216,20 +216,28 @@ GeoHashFix52Bits geohashAlign52Bits(const GeoHashBits hash) {
     return bits;
 }
 
+/* Calculate distance using simplified haversine great circle distance formula.
+ * Given longitude diff is 0 the asin(sqrt(a)) on the haversine is asin(sin(abs(u))).
+ * arcsin(sin(x)) equal to x when x ∈[−𝜋/2,𝜋/2]. Given latitude is between [−𝜋/2,𝜋/2]
+ * we can simplify arcsin(sin(x)) to x.
+ */
+double geohashGetLatDistance(double lat1d, double lat2d) {
+    return EARTH_RADIUS_IN_METERS * fabs(deg_rad(lat2d) - deg_rad(lat1d));
+}
+
 /* Calculate distance using haversine great circle distance formula. */
 double geohashGetDistance(double lon1d, double lat1d, double lon2d, double lat2d) {
-    double lat1r, lon1r, lat2r, lon2r, u, v;
-    lat1r = deg_rad(lat1d);
+    double lat1r, lon1r, lat2r, lon2r, u, v, a;
     lon1r = deg_rad(lon1d);
-    lat2r = deg_rad(lat2d);
     lon2r = deg_rad(lon2d);
-    u = sin((lat2r - lat1r) / 2);
     v = sin((lon2r - lon1r) / 2);
-    double a = u * u;
-    /* if v == 0 we can avoid doing expensive math */
-    if (v != 0.0){
-        a += cos(lat1r) * cos(lat2r) * v * v;
-    }
+    /* if v == 0 we can avoid doing expensive math when lons are practically the same */
+    if (v == 0.0)
+        return geohashGetLatDistance(lat1d, lat2d);
+    lat1r = deg_rad(lat1d);
+    lat2r = deg_rad(lat2d);
+    u = sin((lat2r - lat1r) / 2);
+    a = u * u + cos(lat1r) * cos(lat2r) * v * v;
     return 2.0 * EARTH_RADIUS_IN_METERS * asin(sqrt(a));
 }
 
@@ -259,7 +267,7 @@ int geohashGetDistanceIfInRectangle(double width_m, double height_m, double x1,
                                     double x2, double y2, double *distance) {
     /* latitude distance is less expensive to compute than longitude distance
      * so we check first for the latitude condition */
-    double lat_distance = geohashGetDistance(x2, y2, x2, y1);
+    double lat_distance = geohashGetLatDistance(y2, y1);
     if (lat_distance > height_m/2) {
         return 0;
     }