Skip to content

Instantly share code, notes, and snippets.

@ajfriend
Created November 28, 2025 04:19
Show Gist options
  • Select an option

  • Save ajfriend/ba522a38fa46dbc54875a9042157d2e3 to your computer and use it in GitHub Desktop.

Select an option

Save ajfriend/ba522a38fa46dbc54875a9042157d2e3 to your computer and use it in GitHub Desktop.
Failed attempt to speed up the LatLng loop area calculation
typedef struct {
double s;
double c;
double lng;
} LatLngPre;
static inline LatLngPre precompute_latlng(LatLng x) {
double lat = x.lat / 2.0 + M_PI / 4.0;
LatLngPre llp = {.s = sin(lat), .c = cos(lat), .lng = x.lng};
return llp;
}
static inline double cagnoli_pre(LatLngPre x, LatLngPre y) {
double sa = x.s * y.s;
double ca = x.c * y.c;
double d = y.lng - x.lng;
double sd = sin(d);
double cd = cos(d);
return -2.0 * atan2(sa * sd, sa * cd + ca);
}
static inline double vertArea2(LatLng *verts, int numVerts) {
Kahan k = {0.0, 0.0};
LatLngPre x = precompute_latlng(verts[0]);
LatLngPre first = x;
LatLngPre y;
for (int i = 1; i < numVerts; i++) {
y = precompute_latlng(verts[i]);
kadd(&k, cagnoli_pre(x, y));
x = y;
}
kadd(&k, cagnoli_pre(x, first));
if (k.sum < 0.0) {
kadd(&k, 4.0 * M_PI);
}
return kresult(k);
}
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 d = y.lng - x.lng;
double sd = sin(d);
double cd = cos(d);
return -2.0 * atan2(sa * sd, sa * cd + ca);
}
static inline double vertArea(LatLng *verts, int numVerts) {
Kahan k = {0.0, 0.0};
for (int i = 0; i < numVerts; i++) {
int j = (i + 1) % numVerts;
kadd(&k, cagnoli(verts[i], verts[j]));
}
if (k.sum < 0.0) {
kadd(&k, 4.0 * M_PI);
}
return kresult(k);
}
@ajfriend
Copy link
Author

Since we do cagnoli_pre(verts[i], verts[i+1]), the idea was to cache some of the computation for each point to reuse it, since each point is involved in two computations:

static inline LatLngPre precompute_latlng(LatLng x) {
    double lat = x.lat / 2.0 + M_PI / 4.0;
    LatLngPre llp = {.s = sin(lat), .c = cos(lat), .lng = x.lng};
    return llp;
}

However, this only produced--at most--about 2% improvement. Probably because the remaining trig work, especially the atan2 dominate:

    double sa = x.s * y.s;
    double ca = x.c * y.c;

    double d = y.lng - x.lng;
    double sd = sin(d);
    double cd = cos(d);

    return -2.0 * atan2(sa * sd, sa * cd + ca);
}

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment