diff --git a/README.md b/README.md
index 469660f..8f87355 100644
--- a/README.md
+++ b/README.md
@@ -3,6 +3,7 @@
 [![Swift Version](https://img.shields.io/badge/swift-5.1-blue.svg)](https://swift.org) 
 ![Platform](https://img.shields.io/badge/platform-osx--64|linux--64-lightgrey.svg)
 [![Build Travis-CI Status](https://travis-ci.org/dastrobu/geodesic.svg?branch=master)](https://travis-ci.org/dastrobu/geodesic) 
+[![Swift Version](https://img.shields.io/badge/GeographicLib-1.50.1-blue.svg)](https://geographiclib.sourceforge.io/) 
 
 Solver for the inverse geodesic problem in Swift.
 
@@ -24,8 +25,8 @@ and that's it.
 ## Implementation Details
 
 This Swift package is a wrapper for the 
-[C library for Geodesics 1.49](https://geographiclib.sourceforge.io/html/C/).
-The author of this library is Charles F. F. Karney (charles@karney.com). 
+[C library for Geodesics](https://geographiclib.sourceforge.io/html/C/).
+The author of this library is Charles Karney (charles@karney.com). 
 The goal of this Swift package is to make some algorithms from 
 [GeographicLib](https://geographiclib.sourceforge.io/) available to the Swift world.
 Alternatively one can employ the package 
diff --git a/Sources/geographiclib/LICENSE b/Sources/geographiclib/LICENSE.txt
similarity index 92%
rename from Sources/geographiclib/LICENSE
rename to Sources/geographiclib/LICENSE.txt
index 485d696..1fca9a5 100644
--- a/Sources/geographiclib/LICENSE
+++ b/Sources/geographiclib/LICENSE.txt
@@ -1,9 +1,7 @@
-This is the original LICENSE from GeographicLib:
-
 The MIT License (MIT); this license applies to GeographicLib,
 versions 1.12 and later.
 
-Copyright (c) 2008-2017, Charles Karney
+Copyright (c) 2008-2019, Charles Karney
 
 Permission is hereby granted, free of charge, to any person
 obtaining a copy of this software and associated documentation
diff --git a/Sources/geographiclib/geodesic.c b/Sources/geographiclib/geodesic.c
index aa5a26d..2be842e 100644
--- a/Sources/geographiclib/geodesic.c
+++ b/Sources/geographiclib/geodesic.c
@@ -18,17 +18,28 @@
  *
  * See the comments in geodesic.h for documentation.
  *
- * Copyright (c) Charles Karney (2012-2017) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2019) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * https://geographiclib.sourceforge.io/
  */
 
 #include "geodesic.h"
 #include <math.h>
+#include <limits.h>
+#include <float.h>
 
 #if !defined(HAVE_C99_MATH)
+#if defined(PROJ_LIB)
+/* PROJ requires C99 so HAVE_C99_MATH is implicit */
+#define HAVE_C99_MATH 1
+#else
 #define HAVE_C99_MATH 0
 #endif
+#endif
+
+#if !defined(__cplusplus)
+#define nullptr 0
+#endif
 
 #define GEOGRAPHICLIB_GEODESIC_ORDER 6
 #define nA1   GEOGRAPHICLIB_GEODESIC_ORDER
@@ -56,21 +67,9 @@ static real epsilon, realmin, pi, degree, NaN,
 
 static void Init() {
     if (!init) {
-#if defined(__DBL_MANT_DIG__)
-        digits = __DBL_MANT_DIG__;
-#else
-        digits = 53;
-#endif
-#if defined(__DBL_EPSILON__)
-        epsilon = __DBL_EPSILON__;
-#else
-        epsilon = pow(0.5, digits - 1);
-#endif
-#if defined(__DBL_MIN__)
-        realmin = __DBL_MIN__;
-#else
-        realmin = pow(0.5, 1022);
-#endif
+        digits = DBL_MANT_DIG;
+        epsilon = DBL_EPSILON;
+        realmin = DBL_MIN;
 #if defined(M_PI)
         pi = M_PI;
 #else
@@ -89,10 +88,15 @@ static void Init() {
         tolb = tol0 * tol2;
         xthresh = 1000 * tol2;
         degree = pi / 180;
+#if defined(NAN)
+        NaN = NAN;                  /* NAN is defined in C99 */
+#else
         {
-            real minus1 = -1;
-            NaN = sqrt(minus1);
+          real minus1 = -1;
+          /* cppcheck-suppress wrongmathcall */
+          NaN = sqrt(minus1);
         }
+#endif
         init = 1;
     }
 }
@@ -108,17 +112,30 @@ enum captype {
     OUT_ALL = 0x7F80U
 };
 
-static real sq(real x) {
-    return x * x;
-}
-
 #if HAVE_C99_MATH
+#define hypotx hypot
+/* no need to redirect log1px, since it's only used by atanhx */
 #define atanhx atanh
 #define copysignx copysign
-#define hypotx hypot
 #define cbrtx cbrt
+#define remainderx remainder
+#define remquox remquo
 #else
 
+/* Replacements for C99 math functions */
+
+static real hypotx(real x, real y) {
+    x = fabs(x);
+    y = fabs(y);
+    if (x < y) {
+        x /= y;                     /* y is nonzero */
+        return y * sqrt(1 + x * x);
+    } else {
+        y /= (x != 0 ? x : 1);
+        return x * sqrt(1 + y * y);
+    }
+}
+
 static real log1px(real x) {
     volatile real
             y = 1 + x,
@@ -133,24 +150,63 @@ static real log1px(real x) {
 static real atanhx(real x) {
     real y = fabs(x);             /* Enforce odd parity */
     y = log1px(2 * y / (1 - y)) / 2;
-    return x < 0 ? -y : y;
+    return x > 0 ? y : (x < 0 ? -y : x); /* atanh(-0.0) = -0.0 */
 }
 
 static real copysignx(real x, real y) {
+    /* 1/y trick to get the sign of -0.0 */
     return fabs(x) * (y < 0 || (y == 0 && 1 / y < 0) ? -1 : 1);
 }
 
-static real hypotx(real x, real y) {
-    return sqrt(x * x + y * y);
-}
-
 static real cbrtx(real x) {
-    real y = pow(fabs(x), 1 / (real) (3)); /* Return the real cube root */
-    return x < 0 ? -y : y;
+    real y = pow(fabs(x), 1 / (real) (3));  /* Return the real cube root */
+    return x > 0 ? y : (x < 0 ? -y : x); /* cbrt(-0.0) = -0.0 */
+}
+
+static real remainderx(real x, real y) {
+    real z;
+    y = fabs(y);                 /* The result doesn't depend on the sign of y */
+    z = fmod(x, y);
+    if (z == 0)
+        /* This shouldn't be necessary.  However, before version 14 (2015),
+         * Visual Studio had problems dealing with -0.0.  Specifically
+         *   VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
+         * python 2.7 on Windows 32-bit machines has the same problem. */
+        z = copysignx(z, x);
+    else if (2 * fabs(z) == y)
+        z -= fmod(x, 2 * y) - z;    /* Implement ties to even */
+    else if (2 * fabs(z) > y)
+        z += (z < 0 ? y : -y);      /* Fold remaining cases to (-y/2, y/2) */
+    return z;
+}
+
+static real remquox(real x, real y, int *n) {
+    real z = remainderx(x, y);
+    if (n) {
+        real
+                a = remainderx(x, 2 * y),
+                b = remainderx(x, 4 * y),
+                c = remainderx(x, 8 * y);
+        *n = (a > z ? 1 : (a < z ? -1 : 0));
+        *n += (b > a ? 2 : (b < a ? -2 : 0));
+        *n += (c > b ? 4 : (c < b ? -4 : 0));
+        if (y < 0) *n *= -1;
+        if (y != 0) {
+            if (x / y > 0 && *n <= 0)
+                *n += 8;
+            else if (x / y < 0 && *n >= 0)
+                *n -= 8;
+        }
+    }
+    return z;
 }
 
 #endif
 
+static real sq(real x) {
+    return x * x;
+}
+
 static real sumx(real u, real v, real *t) {
     volatile real s = u + v;
     volatile real up = s - v;
@@ -192,13 +248,8 @@ static void norm2(real *sinx, real *cosx) {
 }
 
 static real AngNormalize(real x) {
-#if HAVE_C99_MATH
-    x = remainder(x, (real)(360));
+    x = remainderx(x, (real) (360));
     return x != -180 ? x : 180;
-#else
-    x = fmod(x, (real) (360));
-    return x <= -180 ? x + 360 : (x <= 180 ? x : x - 360);
-#endif
 }
 
 static real LatFix(real x) {
@@ -231,20 +282,23 @@ static void sincosdx(real x, real *sinx, real *cosx) {
      * the argument to the range [-45, 45] before converting it to radians. */
     real r, s, c;
     int q;
-#if HAVE_C99_MATH && !defined(__GNUC__)
-    /* Disable for gcc because of bug in glibc version < 2.22, see
-     * https://sourceware.org/bugzilla/show_bug.cgi?id=17569 */
-    r = remquo(x, (real)(90), &q);
-#else
-    r = fmod(x, (real) (360));
-    q = (int) (floor(r / 90 + (real) (0.5)));
-    r -= 90 * q;
-#endif
+    r = remquox(x, (real) (90), &q);
     /* now abs(r) <= 45 */
     r *= degree;
     /* Possibly could call the gnu extension sincos */
     s = sin(r);
     c = cos(r);
+#if defined(_MSC_VER) && _MSC_VER < 1900
+    /*
+     * Before version 14 (2015), Visual Studio had problems dealing
+     * with -0.0.  Specifically
+     *   VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
+     *   VC 12       and 64-bit compile:  sin(-0.0)        -> +0.0
+     * AngNormalize has a similar fix.
+     * python 2.7 on Windows 32-bit machines has the same problem.
+     */
+    if (x == 0) s = x;
+#endif
     switch ((unsigned) q & 3U) {
         case 0U:
             *sinx = s;
@@ -349,7 +403,7 @@ static real Lambda12(const struct geod_geodesic *g,
         real *pssig1, real *pcsig1,
         real *pssig2, real *pcsig2,
         real *peps,
-        real *pgomg12,
+        real *pdomg12,
         boolx diffp, real *pdlam12,
         /* Scratch area of the right size */
         real Ca[]);
@@ -384,6 +438,14 @@ static real accsum(const real s[], real y);
 
 static void accneg(real s[]);
 
+static void accrem(real s[], real y);
+
+static real areareduceA(real area[], real area0,
+        int crossings, boolx reverse, boolx sign);
+
+static real areareduceB(real area, real area0,
+        int crossings, boolx reverse, boolx sign);
+
 void geod_init(struct geod_geodesic *g, real a, real f) {
     if (!init) Init();
     g->a = a;
@@ -518,10 +580,10 @@ void geod_lineinit(struct geod_geodesicline *l,
 void geod_gendirectline(struct geod_geodesicline *l,
         const struct geod_geodesic *g,
         real lat1, real lon1, real azi1,
-        unsigned flags, real a12_s12,
+        unsigned flags, real s12_a12,
         unsigned caps) {
     geod_lineinit(l, g, lat1, lon1, azi1, caps);
-    geod_gensetdistance(l, flags, a12_s12);
+    geod_gensetdistance(l, flags, s12_a12);
 }
 
 void geod_directline(struct geod_geodesicline *l,
@@ -544,13 +606,13 @@ real geod_genposition(const struct geod_geodesicline *l,
     real omg12, lam12, lon12;
     real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2, dn2;
     unsigned outmask =
-            (plat2 ? GEOD_LATITUDE : 0U) |
-                    (plon2 ? GEOD_LONGITUDE : 0U) |
-                    (pazi2 ? GEOD_AZIMUTH : 0U) |
-                    (ps12 ? GEOD_DISTANCE : 0U) |
-                    (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
-                    (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
-                    (pS12 ? GEOD_AREA : 0U);
+            (plat2 ? GEOD_LATITUDE : GEOD_NONE) |
+                    (plon2 ? GEOD_LONGITUDE : GEOD_NONE) |
+                    (pazi2 ? GEOD_AZIMUTH : GEOD_NONE) |
+                    (ps12 ? GEOD_DISTANCE : GEOD_NONE) |
+                    (pm12 ? GEOD_REDUCEDLENGTH : GEOD_NONE) |
+                    (pM12 || pM21 ? GEOD_GEODESICSCALE : GEOD_NONE) |
+                    (pS12 ? GEOD_AREA : GEOD_NONE);
 
     outmask &= l->caps & OUT_ALL;
     if (!(TRUE /*Init()*/ &&
@@ -631,7 +693,7 @@ real geod_genposition(const struct geod_geodesicline *l,
     calp2 = l->calp0 * csig2; /* No need to normalize */
 
     if (outmask & GEOD_DISTANCE)
-        s12 = flags & GEOD_ARCMODE ?
+        s12 = (flags & GEOD_ARCMODE) ?
                 l->b * ((1 + l->A1m1) * sig12 + AB1) :
                 s12_a12;
 
@@ -641,7 +703,7 @@ real geod_genposition(const struct geod_geodesicline *l,
         somg2 = l->salp0 * ssig2;
         comg2 = csig2;  /* No need to normalize */
         /* omg12 = omg2 - omg1 */
-        omg12 = flags & GEOD_LONG_UNROLL
+        omg12 = (flags & GEOD_LONG_UNROLL)
                 ? E * (sig12
                         - (atan2(ssig2, csig2) - atan2(l->ssig1, l->csig1))
                         + (atan2(E * somg2, comg2) - atan2(E * l->somg1, l->comg1)))
@@ -651,7 +713,7 @@ real geod_genposition(const struct geod_geodesicline *l,
                 (sig12 + (SinCosSeries(TRUE, ssig2, csig2, l->C3a, nC3 - 1)
                         - l->B31));
         lon12 = lam12 / degree;
-        lon2 = flags & GEOD_LONG_UNROLL ? l->lon1 + lon12 :
+        lon2 = (flags & GEOD_LONG_UNROLL) ? l->lon1 + lon12 :
                 AngNormalize(AngNormalize(l->lon1) + AngNormalize(lon12));
     }
 
@@ -704,47 +766,58 @@ real geod_genposition(const struct geod_geodesicline *l,
         S12 = l->c2 * atan2(salp12, calp12) + l->A4 * (B42 - l->B41);
     }
 
-    if (outmask & GEOD_LATITUDE)
+    /* In the pattern
+     *
+     *   if ((outmask & GEOD_XX) && pYY)
+     *     *pYY = YY;
+     *
+     * the second check "&& pYY" is redundant.  It's there to make the CLang
+     * static analyzer happy.
+     */
+    if ((outmask & GEOD_LATITUDE) && plat2)
         *plat2 = lat2;
-    if (outmask & GEOD_LONGITUDE)
+    if ((outmask & GEOD_LONGITUDE) && plon2)
         *plon2 = lon2;
-    if (outmask & GEOD_AZIMUTH)
+    if ((outmask & GEOD_AZIMUTH) && pazi2)
         *pazi2 = azi2;
-    if (outmask & GEOD_DISTANCE)
+    if ((outmask & GEOD_DISTANCE) && ps12)
         *ps12 = s12;
-    if (outmask & GEOD_REDUCEDLENGTH)
+    if ((outmask & GEOD_REDUCEDLENGTH) && pm12)
         *pm12 = m12;
     if (outmask & GEOD_GEODESICSCALE) {
         if (pM12) *pM12 = M12;
         if (pM21) *pM21 = M21;
     }
-    if (outmask & GEOD_AREA)
+    if ((outmask & GEOD_AREA) && pS12)
         *pS12 = S12;
 
-    return flags & GEOD_ARCMODE ? s12_a12 : sig12 / degree;
+    return (flags & GEOD_ARCMODE) ? s12_a12 : sig12 / degree;
 }
 
 void geod_setdistance(struct geod_geodesicline *l, real s13) {
     l->s13 = s13;
-    l->a13 = geod_genposition(l, GEOD_NOFLAGS, l->s13, 0, 0, 0, 0, 0, 0, 0, 0);
+    l->a13 = geod_genposition(l, GEOD_NOFLAGS, l->s13, nullptr, nullptr, nullptr,
+            nullptr, nullptr, nullptr, nullptr, nullptr);
 }
 
 static void geod_setarc(struct geod_geodesicline *l, real a13) {
     l->a13 = a13;
     l->s13 = NaN;
-    geod_genposition(l, GEOD_ARCMODE, l->a13, 0, 0, 0, &l->s13, 0, 0, 0, 0);
+    geod_genposition(l, GEOD_ARCMODE, l->a13, nullptr, nullptr, nullptr, &l->s13,
+            nullptr, nullptr, nullptr, nullptr);
 }
 
 void geod_gensetdistance(struct geod_geodesicline *l,
         unsigned flags, real s13_a13) {
-    flags & GEOD_ARCMODE ?
+    (flags & GEOD_ARCMODE) ?
             geod_setarc(l, s13_a13) :
             geod_setdistance(l, s13_a13);
 }
 
 void geod_position(const struct geod_geodesicline *l, real s12,
         real *plat2, real *plon2, real *pazi2) {
-    geod_genposition(l, FALSE, s12, plat2, plon2, pazi2, 0, 0, 0, 0, 0);
+    geod_genposition(l, FALSE, s12, plat2, plon2, pazi2,
+            nullptr, nullptr, nullptr, nullptr, nullptr);
 }
 
 real geod_gendirect(const struct geod_geodesic *g,
@@ -755,18 +828,18 @@ real geod_gendirect(const struct geod_geodesic *g,
         real *pS12) {
     struct geod_geodesicline l;
     unsigned outmask =
-            (plat2 ? GEOD_LATITUDE : 0U) |
-                    (plon2 ? GEOD_LONGITUDE : 0U) |
-                    (pazi2 ? GEOD_AZIMUTH : 0U) |
-                    (ps12 ? GEOD_DISTANCE : 0U) |
-                    (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
-                    (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
-                    (pS12 ? GEOD_AREA : 0U);
+            (plat2 ? GEOD_LATITUDE : GEOD_NONE) |
+                    (plon2 ? GEOD_LONGITUDE : GEOD_NONE) |
+                    (pazi2 ? GEOD_AZIMUTH : GEOD_NONE) |
+                    (ps12 ? GEOD_DISTANCE : GEOD_NONE) |
+                    (pm12 ? GEOD_REDUCEDLENGTH : GEOD_NONE) |
+                    (pM12 || pM21 ? GEOD_GEODESICSCALE : GEOD_NONE) |
+                    (pS12 ? GEOD_AREA : GEOD_NONE);
 
     geod_lineinit(&l, g, lat1, lon1, azi1,
             /* Automatically supply GEOD_DISTANCE_IN if necessary */
             outmask |
-                    (flags & GEOD_ARCMODE ? GEOD_NONE : GEOD_DISTANCE_IN));
+                    ((flags & GEOD_ARCMODE) ? GEOD_NONE : GEOD_DISTANCE_IN));
     return geod_genposition(&l, flags, s12_a12,
             plat2, plon2, pazi2, ps12, pm12, pM12, pM21, pS12);
 }
@@ -776,7 +849,7 @@ void geod_direct(const struct geod_geodesic *g,
         real s12,
         real *plat2, real *plon2, real *pazi2) {
     geod_gendirect(g, lat1, lon1, azi1, GEOD_NOFLAGS, s12, plat2, plon2, pazi2,
-            0, 0, 0, 0, 0);
+            nullptr, nullptr, nullptr, nullptr, nullptr);
 }
 
 static real geod_geninverse_int(const struct geod_geodesic *g,
@@ -798,10 +871,10 @@ static real geod_geninverse_int(const struct geod_geodesic *g,
     real omg12 = 0, somg12 = 2, comg12 = 0;
 
     unsigned outmask =
-            (ps12 ? GEOD_DISTANCE : 0U) |
-                    (pm12 ? GEOD_REDUCEDLENGTH : 0U) |
-                    (pM12 || pM21 ? GEOD_GEODESICSCALE : 0U) |
-                    (pS12 ? GEOD_AREA : 0U);
+            (ps12 ? GEOD_DISTANCE : GEOD_NONE) |
+                    (pm12 ? GEOD_REDUCEDLENGTH : GEOD_NONE) |
+                    (pM12 || pM21 ? GEOD_GEODESICSCALE : GEOD_NONE) |
+                    (pS12 ? GEOD_AREA : GEOD_NONE);
 
     outmask &= OUT_ALL;
     /* Compute longitude difference (AngDiff does this carefully).  Result is
@@ -900,9 +973,9 @@ static real geod_geninverse_int(const struct geod_geodesic *g,
         sig12 = atan2(maxx((real) (0), csig1 * ssig2 - ssig1 * csig2),
                 csig1 * csig2 + ssig1 * ssig2);
         Lengths(g, g->n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
-                cbet1, cbet2, &s12x, &m12x, 0,
-                outmask & GEOD_GEODESICSCALE ? &M12 : 0,
-                outmask & GEOD_GEODESICSCALE ? &M21 : 0,
+                cbet1, cbet2, &s12x, &m12x, nullptr,
+                (outmask & GEOD_GEODESICSCALE) ? &M12 : nullptr,
+                (outmask & GEOD_GEODESICSCALE) ? &M21 : nullptr,
                 Ca);
         /* Add the check for sig12 since zero length geodesics might yield m12 <
          * 0.  Test case was
@@ -975,8 +1048,9 @@ static real geod_geninverse_int(const struct geod_geodesic *g,
             unsigned numit = 0;
             /* Bracketing range */
             real salp1a = tiny, calp1a = 1, salp1b = tiny, calp1b = -1;
-            boolx tripn, tripb;
-            for (tripn = FALSE, tripb = FALSE; numit < maxit2; ++numit) {
+            boolx tripn = FALSE;
+            boolx tripb = FALSE;
+            for (; numit < maxit2; ++numit) {
                 /* the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
                  * WGS84 and random input: mean = 2.85, sd = 0.60 */
                 real dv = 0,
@@ -1029,9 +1103,9 @@ static real geod_geninverse_int(const struct geod_geodesic *g,
                         fabs(salp1 - salp1b) + (calp1 - calp1b) < tolb);
             }
             Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
-                    cbet1, cbet2, &s12x, &m12x, 0,
-                    outmask & GEOD_GEODESICSCALE ? &M12 : 0,
-                    outmask & GEOD_GEODESICSCALE ? &M21 : 0, Ca);
+                    cbet1, cbet2, &s12x, &m12x, nullptr,
+                    (outmask & GEOD_GEODESICSCALE) ? &M12 : nullptr,
+                    (outmask & GEOD_GEODESICSCALE) ? &M21 : nullptr, Ca);
             m12x *= g->b;
             s12x *= g->b;
             a12 = sig12 / degree;
@@ -1164,9 +1238,9 @@ void geod_inverseline(struct geod_geodesicline *l,
         real lat1, real lon1, real lat2, real lon2,
         unsigned caps) {
     real salp1, calp1,
-            a12 = geod_geninverse_int(g, lat1, lon1, lat2, lon2, 0,
-            &salp1, &calp1, 0, 0,
-            0, 0, 0, 0),
+            a12 = geod_geninverse_int(g, lat1, lon1, lat2, lon2, nullptr,
+            &salp1, &calp1, nullptr, nullptr,
+            nullptr, nullptr, nullptr, nullptr),
             azi1 = atan2dx(salp1, calp1);
     caps = caps ? caps : GEOD_DISTANCE_IN | GEOD_LONGITUDE;
     /* Ensure that a12 can be converted to a distance */
@@ -1178,7 +1252,8 @@ void geod_inverseline(struct geod_geodesicline *l,
 void geod_inverse(const struct geod_geodesic *g,
         real lat1, real lon1, real lat2, real lon2,
         real *ps12, real *pazi1, real *pazi2) {
-    geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2, 0, 0, 0, 0);
+    geod_geninverse(g, lat1, lon1, lat2, lon2, ps12, pazi1, pazi2,
+            nullptr, nullptr, nullptr, nullptr);
 }
 
 real SinCosSeries(boolx sinp, real sinx, real cosx, const real c[], int n) {
@@ -1190,8 +1265,8 @@ real SinCosSeries(boolx sinp, real sinx, real cosx, const real c[], int n) {
     real ar, y0, y1;
     c += (n + sinp);              /* Point to one beyond last element */
     ar = 2 * (cosx - sinx) * (cosx + sinx); /* 2 * cos(2 * x) */
-    y0 = n & 1 ? *--c : 0;
-    y1 = 0;          /* accumulators for sum */
+    y0 = (n & 1) ? *--c : 0;
+    y1 = 0;        /* accumulators for sum */
     /* Now n is even */
     n /= 2;
     while (n--) {
@@ -1342,22 +1417,7 @@ real InverseStart(const struct geod_geodesic *g,
     boolx shortline = cbet12 >= 0 && sbet12 < (real) (0.5) &&
             cbet2 * lam12 < (real) (0.5);
     real somg12, comg12, ssig12, csig12;
-#if defined(__GNUC__) && __GNUC__ == 4 && \
-  (__GNUC_MINOR__ < 6 || defined(__MINGW32__))
-    /* Volatile declaration needed to fix inverse cases
-     * 88.202499451857 0 -88.202499451857 179.981022032992859592
-     * 89.262080389218 0 -89.262080389218 179.992207982775375662
-     * 89.333123580033 0 -89.333123580032997687 179.99295812360148422
-     * which otherwise fail with g++ 4.4.4 x86 -O3 (Linux)
-     * and g++ 4.4.0 (mingw) and g++ 4.6.1 (tdm mingw). */
-    {
-        volatile real xx1 = sbet2 * cbet1;
-        volatile real xx2 = cbet2 * sbet1;
-        sbet12a = xx1 + xx2;
-    }
-#else
     sbet12a = sbet2 * cbet1 + cbet2 * sbet1;
-#endif
     if (shortline) {
         real sbetm2 = sq(sbet1 + sbet2), omg12;
         /* sin((bet1+bet2)/2)^2
@@ -1423,7 +1483,7 @@ real InverseStart(const struct geod_geodesic *g,
              * Inverse. */
             Lengths(g, g->n, pi + bet12a,
                     sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
-                    cbet1, cbet2, 0, &m12b, &m0, 0, 0, Ca);
+                    cbet1, cbet2, nullptr, &m12b, &m0, nullptr, nullptr, Ca);
             x = -1 + m12b / (cbet1 * cbet2 * m0 * pi);
             betscale = x < -(real) (0.01) ? sbet12a / x :
                     -g->f * sq(cbet1) * pi;
@@ -1588,7 +1648,7 @@ real Lambda12(const struct geod_geodesic *g,
             dlam12 = -2 * g->f1 * dn1 / sbet1;
         else {
             Lengths(g, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2,
-                    cbet1, cbet2, 0, &dlam12, 0, 0, 0, Ca);
+                    cbet1, cbet2, nullptr, &dlam12, nullptr, nullptr, nullptr, Ca);
             dlam12 *= g->f1 / (calp2 * cbet2);
         }
     }
@@ -1876,23 +1936,18 @@ int transit(real lon1, real lon2) {
     /* Compute lon12 the same way as Geodesic::Inverse. */
     lon1 = AngNormalize(lon1);
     lon2 = AngNormalize(lon2);
-    lon12 = AngDiff(lon1, lon2, 0);
+    lon12 = AngDiff(lon1, lon2, nullptr);
     return lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 :
             (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0);
 }
 
 int transitdirect(real lon1, real lon2) {
-#if HAVE_C99_MATH
-    lon1 = remainder(lon1, (real)(720));
-    lon2 = remainder(lon2, (real)(720));
-    return ( (lon2 >= 0 && lon2 < 360 ? 0 : 1) -
-             (lon1 >= 0 && lon1 < 360 ? 0 : 1) );
-#else
-    lon1 = fmod(lon1, (real) (720));
-    lon2 = fmod(lon2, (real) (720));
-    return (((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
-            ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1));
-#endif
+    /* Compute exactly the parity of
+       int(ceil(lon2 / 360)) - int(ceil(lon1 / 360)) */
+    lon1 = remainderx(lon1, (real) (720));
+    lon2 = remainderx(lon2, (real) (720));
+    return ((lon2 <= 0 && lon2 > -360 ? 1 : 0) -
+            (lon1 <= 0 && lon1 > -360 ? 1 : 0));
 }
 
 void accini(real s[]) {
@@ -1930,6 +1985,12 @@ void accneg(real s[]) {
     s[1] = -s[1];
 }
 
+void accrem(real s[], real y) {
+    /* Reduce to [-y/2, y/2]. */
+    s[0] = remainderx(s[0], y);
+    accadd(s, (real) (0));
+}
+
 void geod_polygon_init(struct geod_polygon *p, boolx polylinep) {
     p->polyline = (polylinep != 0);
     geod_polygon_clear(p);
@@ -1952,7 +2013,8 @@ void geod_polygon_addpoint(const struct geod_geodesic *g,
     } else {
         real s12, S12 = 0;       /* Initialize S12 to stop Visual Studio warning */
         geod_geninverse(g, p->lat, p->lon, lat, lon,
-                &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
+                &s12, nullptr, nullptr, nullptr, nullptr, nullptr,
+                p->polyline ? nullptr : &S12);
         accadd(p->P, s12);
         if (!p->polyline) {
             accadd(p->A, S12);
@@ -1967,11 +2029,14 @@ void geod_polygon_addpoint(const struct geod_geodesic *g,
 void geod_polygon_addedge(const struct geod_geodesic *g,
         struct geod_polygon *p,
         real azi, real s) {
-    if (p->num) {                 /* Do nothing is num is zero */
-        real lat, lon, S12 = 0;  /* Initialize S12 to stop Visual Studio warning */
+    if (p->num) {              /* Do nothing is num is zero */
+        /* Initialize S12 to stop Visual Studio warning.  Initialization of lat and
+         * lon is to make CLang static analyzer happy. */
+        real lat = 0, lon = 0, S12 = 0;
         geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
-                &lat, &lon, 0,
-                0, 0, 0, 0, p->polyline ? 0 : &S12);
+                &lat, &lon, nullptr,
+                nullptr, nullptr, nullptr, nullptr,
+                p->polyline ? nullptr : &S12);
         accadd(p->P, s);
         if (!p->polyline) {
             accadd(p->A, S12);
@@ -1987,8 +2052,7 @@ unsigned geod_polygon_compute(const struct geod_geodesic *g,
         const struct geod_polygon *p,
         boolx reverse, boolx sign,
         real *pA, real *pP) {
-    real s12, S12, t[2], area0;
-    int crossings;
+    real s12, S12, t[2];
     if (p->num < 2) {
         if (pP) *pP = 0;
         if (!p->polyline && pA) *pA = 0;
@@ -1999,31 +2063,14 @@ unsigned geod_polygon_compute(const struct geod_geodesic *g,
         return p->num;
     }
     geod_geninverse(g, p->lat, p->lon, p->lat0, p->lon0,
-            &s12, 0, 0, 0, 0, 0, &S12);
+            &s12, nullptr, nullptr, nullptr, nullptr, nullptr, &S12);
     if (pP) *pP = accsum(p->P, s12);
     acccopy(p->A, t);
     accadd(t, S12);
-    crossings = p->crossings + transit(p->lon, p->lon0);
-    area0 = 4 * pi * g->c2;
-    if (crossings & 1)
-        accadd(t, (t[0] < 0 ? 1 : -1) * area0 / 2);
-    /* area is with the clockwise sense.  If !reverse convert to
-     * counter-clockwise convention. */
-    if (!reverse)
-        accneg(t);
-    /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
-    if (sign) {
-        if (t[0] > area0 / 2)
-            accadd(t, -area0);
-        else if (t[0] <= -area0 / 2)
-            accadd(t, +area0);
-    } else {
-        if (t[0] >= area0)
-            accadd(t, -area0);
-        else if (t[0] < 0)
-            accadd(t, +area0);
-    }
-    if (pA) *pA = 0 + t[0];
+    if (pA)
+        *pA = areareduceA(t, 4 * pi * g->c2,
+                p->crossings + transit(p->lon, p->lon0),
+                reverse, sign);
     return p->num;
 }
 
@@ -2032,7 +2079,7 @@ unsigned geod_polygon_testpoint(const struct geod_geodesic *g,
         real lat, real lon,
         boolx reverse, boolx sign,
         real *pA, real *pP) {
-    real perimeter, tempsum, area0;
+    real perimeter, tempsum;
     int crossings, i;
     unsigned num = p->num + 1;
     if (num == 1) {
@@ -2048,7 +2095,8 @@ unsigned geod_polygon_testpoint(const struct geod_geodesic *g,
         geod_geninverse(g,
                 i == 0 ? p->lat : lat, i == 0 ? p->lon : lon,
                 i != 0 ? p->lat0 : lat, i != 0 ? p->lon0 : lon,
-                &s12, 0, 0, 0, 0, 0, p->polyline ? 0 : &S12);
+                &s12, nullptr, nullptr, nullptr, nullptr, nullptr,
+                p->polyline ? nullptr : &S12);
         perimeter += s12;
         if (!p->polyline) {
             tempsum += S12;
@@ -2061,26 +2109,7 @@ unsigned geod_polygon_testpoint(const struct geod_geodesic *g,
     if (p->polyline)
         return num;
 
-    area0 = 4 * pi * g->c2;
-    if (crossings & 1)
-        tempsum += (tempsum < 0 ? 1 : -1) * area0 / 2;
-    /* area is with the clockwise sense.  If !reverse convert to
-     * counter-clockwise convention. */
-    if (!reverse)
-        tempsum *= -1;
-    /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
-    if (sign) {
-        if (tempsum > area0 / 2)
-            tempsum -= area0;
-        else if (tempsum <= -area0 / 2)
-            tempsum += area0;
-    } else {
-        if (tempsum >= area0)
-            tempsum -= area0;
-        else if (tempsum < 0)
-            tempsum += area0;
-    }
-    if (pA) *pA = 0 + tempsum;
+    if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign);
     return num;
 }
 
@@ -2089,7 +2118,7 @@ unsigned geod_polygon_testedge(const struct geod_geodesic *g,
         real azi, real s,
         boolx reverse, boolx sign,
         real *pA, real *pP) {
-    real perimeter, tempsum, area0;
+    real perimeter, tempsum;
     int crossings;
     unsigned num = p->num + 1;
     if (num == 1) {               /* we don't have a starting point! */
@@ -2106,40 +2135,23 @@ unsigned geod_polygon_testedge(const struct geod_geodesic *g,
     tempsum = p->A[0];
     crossings = p->crossings;
     {
-        real lat, lon, s12, S12;
+        /* Initialization of lat, lon, and S12 is to make CLang static analyzer
+         * happy. */
+        real lat = 0, lon = 0, s12, S12 = 0;
         geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s,
-                &lat, &lon, 0,
-                0, 0, 0, 0, &S12);
+                &lat, &lon, nullptr,
+                nullptr, nullptr, nullptr, nullptr, &S12);
         tempsum += S12;
         crossings += transitdirect(p->lon, lon);
         geod_geninverse(g, lat, lon, p->lat0, p->lon0,
-                &s12, 0, 0, 0, 0, 0, &S12);
+                &s12, nullptr, nullptr, nullptr, nullptr, nullptr, &S12);
         perimeter += s12;
         tempsum += S12;
         crossings += transit(lon, p->lon0);
     }
 
-    area0 = 4 * pi * g->c2;
-    if (crossings & 1)
-        tempsum += (tempsum < 0 ? 1 : -1) * area0 / 2;
-    /* area is with the clockwise sense.  If !reverse convert to
-     * counter-clockwise convention. */
-    if (!reverse)
-        tempsum *= -1;
-    /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
-    if (sign) {
-        if (tempsum > area0 / 2)
-            tempsum -= area0;
-        else if (tempsum <= -area0 / 2)
-            tempsum += area0;
-    } else {
-        if (tempsum >= area0)
-            tempsum -= area0;
-        else if (tempsum < 0)
-            tempsum += area0;
-    }
     if (pP) *pP = perimeter;
-    if (pA) *pA = 0 + tempsum;
+    if (pA) *pA = areareduceB(tempsum, 4 * pi * g->c2, crossings, reverse, sign);
     return num;
 }
 
@@ -2154,4 +2166,52 @@ void geod_polygonarea(const struct geod_geodesic *g,
     geod_polygon_compute(g, &p, FALSE, TRUE, pA, pP);
 }
 
+real areareduceA(real area[], real area0,
+        int crossings, boolx reverse, boolx sign) {
+    accrem(area, area0);
+    if (crossings & 1)
+        accadd(area, (area[0] < 0 ? 1 : -1) * area0 / 2);
+    /* area is with the clockwise sense.  If !reverse convert to
+     * counter-clockwise convention. */
+    if (!reverse)
+        accneg(area);
+    /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+    if (sign) {
+        if (area[0] > area0 / 2)
+            accadd(area, -area0);
+        else if (area[0] <= -area0 / 2)
+            accadd(area, +area0);
+    } else {
+        if (area[0] >= area0)
+            accadd(area, -area0);
+        else if (area[0] < 0)
+            accadd(area, +area0);
+    }
+    return 0 + area[0];
+}
+
+real areareduceB(real area, real area0,
+        int crossings, boolx reverse, boolx sign) {
+    area = remainderx(area, area0);
+    if (crossings & 1)
+        area += (area < 0 ? 1 : -1) * area0 / 2;
+    /* area is with the clockwise sense.  If !reverse convert to
+     * counter-clockwise convention. */
+    if (!reverse)
+        area *= -1;
+    /* If sign put area in (-area0/2, area0/2], else put area in [0, area0) */
+    if (sign) {
+        if (area > area0 / 2)
+            area -= area0;
+        else if (area <= -area0 / 2)
+            area += area0;
+    } else {
+        if (area >= area0)
+            area -= area0;
+        else if (area < 0)
+            area += area0;
+    }
+    return 0 + area;
+}
+
 /** @endcond */
diff --git a/Sources/geographiclib/include/geodesic.h b/Sources/geographiclib/include/geodesic.h
index 8392a6b..731a143 100644
--- a/Sources/geographiclib/include/geodesic.h
+++ b/Sources/geographiclib/include/geodesic.h
@@ -107,12 +107,12 @@
  * twice about restructuring the internals of the C code since this may make
  * porting fixes from the C++ code more difficult.
  *
- * Copyright (c) Charles Karney (2012-2017) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2019) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * https://geographiclib.sourceforge.io/
  *
  * This library was distributed with
- * <a href="../index.html">GeographicLib</a> 1.49.
+ * <a href="../index.html">GeographicLib</a> 1.50.
  **********************************************************************/
 
 #if !defined(GEODESIC_H)
@@ -127,7 +127,7 @@
  * The minor version of the geodesic library.  (This tracks the version of
  * GeographicLib.)
  **********************************************************************/
-#define GEODESIC_VERSION_MINOR 49
+#define GEODESIC_VERSION_MINOR 50
 /**
  * The patch level of the geodesic library.  (This tracks the version of
  * GeographicLib.)
@@ -157,6 +157,20 @@
                       GEODESIC_VERSION_MINOR, \
                       GEODESIC_VERSION_PATCH)
 
+#if !defined(GEOD_DLL)
+#if defined(_MSC_VER)
+#define GEOD_DLL __declspec(dllexport)
+#elif defined(__GNUC__)
+#define GEOD_DLL __attribute__ ((visibility("default")))
+#else
+#define GEOD_DLL
+#endif
+#endif
+
+#if defined(PROJ_RENAME_SYMBOLS)
+#include "proj_symbol_rename.h"
+#endif
+
 #if defined(__cplusplus)
 extern "C" {
 #endif
@@ -224,7 +238,7 @@ struct geod_polygon {
  * @param[in] a the equatorial radius (meters).
  * @param[in] f the flattening.
  **********************************************************************/
-void geod_init(struct geod_geodesic *g, double a, double f);
+void GEOD_DLL geod_init(struct geod_geodesic *g, double a, double f);
 
 /**
  * Solve the direct geodesic problem.
@@ -262,7 +276,7 @@ void geod_init(struct geod_geodesic *g, double a, double f);
  printf("%.5f %.5f\n", lat, lon);
  @endcode
  **********************************************************************/
-void geod_direct(const struct geod_geodesic *g,
+void GEOD_DLL geod_direct(const struct geod_geodesic *g,
         double lat1, double lon1, double azi1, double s12,
         double *plat2, double *plon2, double *pazi2);
 
@@ -304,11 +318,12 @@ void geod_direct(const struct geod_geodesic *g,
  * that the quantity \e lon2 &minus; \e lon1 indicates how many times and in
  * what sense the geodesic encircles the ellipsoid.
  **********************************************************************/
-double geod_gendirect(const struct geod_geodesic *g,
+double GEOD_DLL geod_gendirect(const struct geod_geodesic *g,
         double lat1, double lon1, double azi1,
         unsigned flags, double s12_a12,
         double *plat2, double *plon2, double *pazi2,
-        double *ps12, double *pm12, double *pM12, double *pM21,
+        double *ps12, double *pm12,
+        double *pM12, double *pM21,
         double *pS12);
 
 /**
@@ -349,8 +364,9 @@ double geod_gendirect(const struct geod_geodesic *g,
  printf("%.3f\n", s12);
  @endcode
  **********************************************************************/
-void geod_inverse(const struct geod_geodesic *g,
-        double lat1, double lon1, double lat2, double lon2,
+void GEOD_DLL geod_inverse(const struct geod_geodesic *g,
+        double lat1, double lon1,
+        double lat2, double lon2,
         double *ps12, double *pazi1, double *pazi2);
 
 /**
@@ -380,8 +396,9 @@ void geod_inverse(const struct geod_geodesic *g,
  * "return" arguments \e ps12, etc., may be replaced by 0, if you do not need
  * some quantities computed.
  **********************************************************************/
-double geod_geninverse(const struct geod_geodesic *g,
-        double lat1, double lon1, double lat2, double lon2,
+double GEOD_DLL geod_geninverse(const struct geod_geodesic *g,
+        double lat1, double lon1,
+        double lat2, double lon2,
         double *ps12, double *pazi1, double *pazi2,
         double *pm12, double *pM12, double *pM21,
         double *pS12);
@@ -425,9 +442,10 @@ double geod_geninverse(const struct geod_geodesic *g,
  * When initialized by this function, point 3 is undefined (l->s13 = l->a13 =
  * NaN).
  **********************************************************************/
-void geod_lineinit(struct geod_geodesicline *l,
+void GEOD_DLL geod_lineinit(struct geod_geodesicline *l,
         const struct geod_geodesic *g,
-        double lat1, double lon1, double azi1, unsigned caps);
+        double lat1, double lon1, double azi1,
+        unsigned caps);
 
 /**
  * Initialize a geod_geodesicline object in terms of the direct geodesic
@@ -450,9 +468,10 @@ void geod_lineinit(struct geod_geodesicline *l,
  * 2 of the direct geodesic problem.  See geod_lineinit() for more
  * information.
  **********************************************************************/
-void geod_directline(struct geod_geodesicline *l,
+void GEOD_DLL geod_directline(struct geod_geodesicline *l,
         const struct geod_geodesic *g,
-        double lat1, double lon1, double azi1, double s12,
+        double lat1, double lon1,
+        double azi1, double s12,
         unsigned caps);
 
 /**
@@ -480,7 +499,7 @@ void geod_directline(struct geod_geodesicline *l,
  * 2 of the direct geodesic problem.  See geod_lineinit() for more
  * information.
  **********************************************************************/
-void geod_gendirectline(struct geod_geodesicline *l,
+void GEOD_DLL geod_gendirectline(struct geod_geodesicline *l,
         const struct geod_geodesic *g,
         double lat1, double lon1, double azi1,
         unsigned flags, double s12_a12,
@@ -506,9 +525,10 @@ void geod_gendirectline(struct geod_geodesicline *l,
  * 2 of the inverse geodesic problem.  See geod_lineinit() for more
  * information.
  **********************************************************************/
-void geod_inverseline(struct geod_geodesicline *l,
+void GEOD_DLL geod_inverseline(struct geod_geodesicline *l,
         const struct geod_geodesic *g,
-        double lat1, double lon1, double lat2, double lon2,
+        double lat1, double lon1,
+        double lat2, double lon2,
         unsigned caps);
 
 /**
@@ -524,9 +544,9 @@ void geod_inverseline(struct geod_geodesicline *l,
  * @param[out] pazi2 pointer to the (forward) azimuth at point 2 (degrees).
  *
  * \e l must have been initialized with a call, e.g., to geod_lineinit(),
- * with \e caps |= GEOD_DISTANCE_IN.  The values of \e lon2 and \e azi2
- * returned are in the range [&minus;180&deg;, 180&deg;].  Any of the
- * "return" arguments \e plat2, etc., may be replaced by 0, if you do not
+ * with \e caps |= GEOD_DISTANCE_IN (or \e caps = 0).  The values of \e lon2
+ * and \e azi2 returned are in the range [&minus;180&deg;, 180&deg;].  Any of
+ * the "return" arguments \e plat2, etc., may be replaced by 0, if you do not
  * need some quantities computed.
  *
  * Example, compute way points between JFK and Singapore Changi Airport
@@ -556,7 +576,7 @@ void geod_inverseline(struct geod_geodesicline *l,
  }
  @endcode
  **********************************************************************/
-void geod_position(const struct geod_geodesicline *l, double s12,
+void GEOD_DLL geod_position(const struct geod_geodesicline *l, double s12,
         double *plat2, double *plon2, double *pazi2);
 
 /**
@@ -623,7 +643,7 @@ void geod_position(const struct geod_geodesicline *l, double s12,
  }
  @endcode
  **********************************************************************/
-double geod_genposition(const struct geod_geodesicline *l,
+double GEOD_DLL geod_genposition(const struct geod_geodesicline *l,
         unsigned flags, double s12_a12,
         double *plat2, double *plon2, double *pazi2,
         double *ps12, double *pm12,
@@ -633,19 +653,19 @@ double geod_genposition(const struct geod_geodesicline *l,
 /**
  * Specify position of point 3 in terms of distance.
  *
- * @param[inout] l a pointer to the geod_geodesicline object.
+ * @param[in,out] l a pointer to the geod_geodesicline object.
  * @param[in] s13 the distance from point 1 to point 3 (meters); it
  *   can be negative.
  *
  * This is only useful if the geod_geodesicline object has been constructed
  * with \e caps |= GEOD_DISTANCE_IN.
  **********************************************************************/
-void geod_setdistance(struct geod_geodesicline *l, double s13);
+void GEOD_DLL geod_setdistance(struct geod_geodesicline *l, double s13);
 
 /**
  * Specify position of point 3 in terms of either distance or arc length.
  *
- * @param[inout] l a pointer to the geod_geodesicline object.
+ * @param[in,out] l a pointer to the geod_geodesicline object.
  * @param[in] flags either GEOD_NOFLAGS or GEOD_ARCMODE to determining the
  *   meaning of the \e s13_a13.
  * @param[in] s13_a13 if \e flags = GEOD_NOFLAGS, this is the distance
@@ -657,7 +677,7 @@ void geod_setdistance(struct geod_geodesicline *l, double s13);
  * GEOD_ARCMODE, the \e s13 is only set if the geod_geodesicline object has
  * been constructed with \e caps |= GEOD_DISTANCE.
  **********************************************************************/
-void geod_gensetdistance(struct geod_geodesicline *l,
+void GEOD_DLL geod_gensetdistance(struct geod_geodesicline *l,
         unsigned flags, double s13_a13);
 
 /**
@@ -679,14 +699,14 @@ void geod_gensetdistance(struct geod_geodesicline *l,
  * An example of the use of this function is given in the documentation for
  * geod_polygon_compute().
  **********************************************************************/
-void geod_polygon_init(struct geod_polygon *p, int polylinep);
+void GEOD_DLL geod_polygon_init(struct geod_polygon *p, int polylinep);
 
 /**
  * Clear the polygon, allowing a new polygon to be started.
  *
  * @param[in,out] p a pointer to the object to be cleared.
  **********************************************************************/
-void geod_polygon_clear(struct geod_polygon *p);
+void GEOD_DLL geod_polygon_clear(struct geod_polygon *p);
 
 /**
  * Add a point to the polygon or polyline.
@@ -706,7 +726,7 @@ void geod_polygon_clear(struct geod_polygon *p);
  * An example of the use of this function is given in the documentation for
  * geod_polygon_compute().
  **********************************************************************/
-void geod_polygon_addpoint(const struct geod_geodesic *g,
+void GEOD_DLL geod_polygon_addpoint(const struct geod_geodesic *g,
         struct geod_polygon *p,
         double lat, double lon);
 
@@ -726,7 +746,7 @@ void geod_polygon_addpoint(const struct geod_geodesic *g,
  * added yet.  The \e lat and \e lon fields of \e p give the location of the
  * new vertex.
  **********************************************************************/
-void geod_polygon_addedge(const struct geod_geodesic *g,
+void GEOD_DLL geod_polygon_addedge(const struct geod_geodesic *g,
         struct geod_polygon *p,
         double azi, double s);
 
@@ -749,10 +769,12 @@ void geod_polygon_addedge(const struct geod_geodesic *g,
  *
  * The area and perimeter are accumulated at two times the standard floating
  * point precision to guard against the loss of accuracy with many-sided
- * polygons.  Only simple polygons (which are not self-intersecting) are
- * allowed.  There's no need to "close" the polygon by repeating the first
- * vertex.  Set \e pA or \e pP to zero, if you do not want the corresponding
- * quantity returned.
+ * polygons.  Arbitrarily complex polygons are allowed.  In the case of
+ * self-intersecting polygons the area is accumulated "algebraically", e.g.,
+ * the areas of the 2 loops in a figure-8 polygon will partially cancel.
+ * There's no need to "close" the polygon by repeating the first vertex.  Set
+ * \e pA or \e pP to zero, if you do not want the corresponding quantity
+ * returned.
  *
  * More points can be added to the polygon after this call.
  *
@@ -773,7 +795,7 @@ void geod_polygon_addedge(const struct geod_geodesic *g,
  printf("%d %.8f %.3f\n", n, P, A);
  @endcode
  **********************************************************************/
-unsigned geod_polygon_compute(const struct geod_geodesic *g,
+unsigned GEOD_DLL geod_polygon_compute(const struct geod_geodesic *g,
         const struct geod_polygon *p,
         int reverse, int sign,
         double *pA, double *pP);
@@ -804,7 +826,7 @@ unsigned geod_polygon_compute(const struct geod_geodesic *g,
  *
  * \e lat should be in the range [&minus;90&deg;, 90&deg;].
  **********************************************************************/
-unsigned geod_polygon_testpoint(const struct geod_geodesic *g,
+unsigned GEOD_DLL geod_polygon_testpoint(const struct geod_geodesic *g,
         const struct geod_polygon *p,
         double lat, double lon,
         int reverse, int sign,
@@ -835,7 +857,7 @@ unsigned geod_polygon_testpoint(const struct geod_geodesic *g,
  *   polyline (meters).
  * @return the number of points.
  **********************************************************************/
-unsigned geod_polygon_testedge(const struct geod_geodesic *g,
+unsigned GEOD_DLL geod_polygon_testedge(const struct geod_geodesic *g,
         const struct geod_polygon *p,
         double azi, double s,
         int reverse, int sign,
@@ -854,10 +876,11 @@ unsigned geod_polygon_testedge(const struct geod_geodesic *g,
  *
  * \e lats should be in the range [&minus;90&deg;, 90&deg;].
  *
- * Only simple polygons (which are not self-intersecting) are allowed.
- * There's no need to "close" the polygon by repeating the first vertex.  The
- * area returned is signed with counter-clockwise traversal being treated as
- * positive.
+ * Arbitrarily complex polygons are allowed.  In the case self-intersecting
+ * of polygons the area is accumulated "algebraically", e.g., the areas of
+ * the 2 loops in a figure-8 polygon will partially cancel.  There's no need
+ * to "close" the polygon by repeating the first vertex.  The area returned
+ * is signed with counter-clockwise traversal being treated as positive.
  *
  * Example, compute the area of Antarctica:
  @code{.c}
@@ -873,7 +896,7 @@ unsigned geod_polygon_testedge(const struct geod_geodesic *g,
  printf("%.0f %.2f\n", A, P);
  @endcode
  **********************************************************************/
-void geod_polygonarea(const struct geod_geodesic *g,
+void GEOD_DLL geod_polygonarea(const struct geod_geodesic *g,
         double lats[], double lons[], int n,
         double *pA, double *pP);