Speed up geo coordinate transforms: combined sincos + fewer trig calls

src/h3lib/include/mathExtensions.h

@@ -20,6 +20,7 @@
 #ifndef MATHEXTENSIONS_H
 #define MATHEXTENSIONS_H
 
+#include <math.h>
 #include <stdbool.h>
 #include <stdint.h>
 
@@ -28,6 +29,25 @@
  */
 #define MAX(a, b) (((a) > (b)) ? (a) : (b))
 
+/**
+ * Simultaneously compute the sine and cosine of an angle.
+ *
+ * Many H3 coordinate transforms need both sin(x) and cos(x) of the same angle.
+ * Computing them together lets the math library share the argument reduction
+ * and polynomial evaluation, which is roughly twice as fast as two independent
+ * sin()/cos() calls, while returning identical values. The compiler builtin
+ * lowers to a single libm `sincos` call where one exists and falls back to
+ * separate sin()/cos() otherwise, so this stays portable.
+ */
+static inline void _sincos(double x, double *s, double *c) {
+#if defined(__GNUC__) || defined(__clang__)
+    __builtin_sincos(x, s, c);
+#else
+    *s = sin(x);
+    *c = cos(x);
+#endif
+}
+
 /** Evaluates to true if a + b would overflow for int32 */
 static inline bool ADD_INT32S_OVERFLOWS(int32_t a, int32_t b) {
     if (a > 0) {

src/h3lib/lib/faceijk.c

@@ -30,6 +30,7 @@
 #include "coordijk.h"
 #include "h3Index.h"
 #include "latLng.h"
+#include "mathExtensions.h"
 #include "vec3d.h"
 
 /** square root of 7 and inverse square root of 7 */
@@ -418,9 +419,11 @@ void _geoToHex2d(const LatLng *g, int res, int *face, Vec2d *v) {
 
     // we now have (r, theta) in hex2d with theta ccw from x-axes
 
-    // convert to local x,y
-    v->x = r * cos(theta);
-    v->y = r * sin(theta);
+    // convert to local x,y (both the sine and cosine of theta are needed)
+    double sinTheta, cosTheta;
+    _sincos(theta, &sinTheta, &cosTheta);
+    v->x = r * cosTheta;
+    v->y = r * sinTheta;
 }
 
 /**

src/h3lib/lib/latLng.c

@@ -209,9 +209,15 @@ double H3_EXPORT(greatCircleDistanceM)(const LatLng *a, const LatLng *b) {
  * @return The azimuth in radians from p1 to p2.
  */
 double _geoAzimuthRads(const LatLng *p1, const LatLng *p2) {
-    return atan2(cos(p2->lat) * sin(p2->lng - p1->lng),
-                 cos(p1->lat) * sin(p2->lat) -
-                     sin(p1->lat) * cos(p2->lat) * cos(p2->lng - p1->lng));
+    // Compute each sin/cos pair together and reuse cos(p2->lat) and the
+    // longitude difference instead of evaluating them twice.
+    double sinP1Lat, cosP1Lat, sinP2Lat, cosP2Lat, sinDLng, cosDLng;
+    _sincos(p1->lat, &sinP1Lat, &cosP1Lat);
+    _sincos(p2->lat, &sinP2Lat, &cosP2Lat);
+    _sincos(p2->lng - p1->lng, &sinDLng, &cosDLng);
+
+    return atan2(cosP2Lat * sinDLng,
+                 cosP1Lat * sinP2Lat - sinP1Lat * cosP2Lat * cosDLng);
 }
 
 /**
@@ -254,8 +260,14 @@ void _geoAzDistanceRads(const LatLng *p1, double az, double distance,
             p2->lng = constrainLng(p1->lng);
     } else  // not due north or south
     {
-        sinlat = sin(p1->lat) * cos(distance) +
-                 cos(p1->lat) * sin(distance) * cos(az);
+        // Each of p1->lat, the distance, and the azimuth is used for both its
+        // sine and its cosine, so compute the pairs together and reuse them.
+        double sinP1Lat, cosP1Lat, sinDist, cosDist, sinAz, cosAz;
+        _sincos(p1->lat, &sinP1Lat, &cosP1Lat);
+        _sincos(distance, &sinDist, &cosDist);
+        _sincos(az, &sinAz, &cosAz);
+
+        sinlat = sinP1Lat * cosDist + cosP1Lat * sinDist * cosAz;
         if (sinlat > 1.0) sinlat = 1.0;
         if (sinlat < -1.0) sinlat = -1.0;
         p2->lat = asin(sinlat);
@@ -269,9 +281,9 @@ void _geoAzDistanceRads(const LatLng *p1, double az, double distance,
             p2->lng = 0.0;
         } else {
             double invcosp2lat = 1.0 / cos(p2->lat);
-            sinlng = sin(az) * sin(distance) * invcosp2lat;
-            coslng = (cos(distance) - sin(p1->lat) * sin(p2->lat)) /
-                     cos(p1->lat) * invcosp2lat;
+            sinlng = sinAz * sinDist * invcosp2lat;
+            coslng =
+                (cosDist - sinP1Lat * sin(p2->lat)) / cosP1Lat * invcosp2lat;
             if (sinlng > 1.0) sinlng = 1.0;
             if (sinlng < -1.0) sinlng = -1.0;
             if (coslng > 1.0) coslng = 1.0;

src/h3lib/lib/vec3d.c

@@ -21,6 +21,8 @@
 
 #include <math.h>
 
+#include "mathExtensions.h"
+
 /**
  * Square of a number
  *
@@ -48,9 +50,15 @@ double _pointSquareDist(const Vec3d *v1, const Vec3d *v2) {
  * @param v The 3D coordinate of the point.
  */
 void _geoToVec3d(const LatLng *geo, Vec3d *v) {
-    double r = cos(geo->lat);
+    // sin and cos of the latitude and of the longitude are each needed as a
+    // pair, so compute them together instead of with four separate calls.
+    double sinLat, cosLat, sinLng, cosLng;
+    _sincos(geo->lat, &sinLat, &cosLat);
+    _sincos(geo->lng, &sinLng, &cosLng);
+
+    double r = cosLat;
 
-    v->z = sin(geo->lat);
-    v->x = cos(geo->lng) * r;
-    v->y = sin(geo->lng) * r;
+    v->z = sinLat;
+    v->x = cosLng * r;
+    v->y = sinLng * r;
 }