Make _sincos reach libm's sincos on every Linux libc, and use it in cagnoli

src/h3lib/include/mathExtensions.h

@@ -35,12 +35,26 @@
  * 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.
+ * sin()/cos() calls, while returning identical values.
+ *
+ * On Linux the libm `sincos` is called directly. `__builtin_sincos`, a
+ * `sin`/`cos` pair the optimizer merges by itself, and even a plain call to
+ * a function named `sincos` all become the `llvm.sincos` intrinsic, which
+ * LLVM lowers to a `sincos` call only for the C libraries it knows to have
+ * one (glibc, Android, Fuchsia) and splits back into separate `sin` and
+ * `cos` calls for the others, musl among them. Declaring the function under
+ * another C name with the `sincos` assembler label keeps the compiler from
+ * recognizing it, so the call reaches libm as written. `sincos` is a GNU
+ * extension that every Linux libc provides. Elsewhere the builtin is used.
  */
+#if defined(__linux__)
+void _h3Sincos(double x, double *s, double *c) __asm__("sincos");
+#endif
+
 static inline void _sincos(double x, double *s, double *c) {
-#if defined(__GNUC__) || defined(__clang__)
+#if defined(__linux__)
+    _h3Sincos(x, s, c);
+#elif defined(__GNUC__) || defined(__clang__)
     __builtin_sincos(x, s, c);
 #else
     *s = sin(x);

src/h3lib/lib/area.c

@@ -26,6 +26,7 @@
 #include "constants.h"
 #include "h3Assert.h"
 #include "h3api.h"
+#include "mathExtensions.h"
 
 // Cagnoli contribution for edge arc x to y, following d3-geo’s
 // area implementation:
@@ -34,8 +35,11 @@ static inline double cagnoli(LatLng x, LatLng y) {
     x.lat = x.lat / 2.0 + M_PI / 4.0;
     y.lat = y.lat / 2.0 + M_PI / 4.0;
 
-    double sa = sin(x.lat) * sin(y.lat);
-    double ca = cos(x.lat) * cos(y.lat);
+    double sinXLat, cosXLat, sinYLat, cosYLat;
+    _sincos(x.lat, &sinXLat, &cosXLat);
+    _sincos(y.lat, &sinYLat, &cosYLat);
+    double sa = sinXLat * sinYLat;
+    double ca = cosXLat * cosYLat;
 
     // For an edge crossing the antimeridian, the raw difference is close to
     // +-2*pi, and the true (tiny) longitude step is lost to the ulp of 2*pi.
@@ -46,8 +50,8 @@ static inline double cagnoli(LatLng x, LatLng y) {
     if (d > M_PI || d < -M_PI) {
         d = (y.lng - copysign(M_PI, y.lng)) - (x.lng - copysign(M_PI, x.lng));
     }
-    double sd = sin(d);
-    double cd = cos(d);
+    double sd, cd;
+    _sincos(d, &sd, &cd);
 
     return -2.0 * atan2(sa * sd, sa * cd + ca);
 }