gie.c and builtins.gie now able to reproduce internal test results

improved docs, improved strtod - avoid precision loss for very long fractions

Switch gie.c to use same framework as cct.c

numerous improvements in proj_strtod.c and gie.c

Add gie to the build system

7 files changed, 823 insertions, 30 deletions

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 788273a9..494eef9b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,11 +3,12 @@
 include(lib_proj.cmake)
 
 # configure executable build
-option(BUILD_CCT   "Build cct (coordinate conversion and transformation tool)" ON)
-option(BUILD_CS2CS "Build cs2cs (coordinate systems to coordinate systems translation tool)" ON)
-option(BUILD_PROJ "Build proj (cartographic projection tool : latlong <-> projected coordinates" ON)
-option(BUILD_GEOD "Build geod (computation of geodesic lines)" ON)
-option(BUILD_NAD2BIN "Build nad2bin (format conversion tool) " ON)
+option(BUILD_CCT     "Build cct (coordinate conversion and transformation tool)" ON)
+option(BUILD_CS2CS   "Build cs2cs (coordinate systems to coordinate systems translation tool)" ON)
+option(BUILD_GEOD    "Build geod (computation of geodesic lines)" ON)
+option(BUILD_GIE     "Build gie (geospatial integrity investigation environment - a PROJ.4 test tool)" ON)
+option(BUILD_NAD2BIN "Build nad2bin (format conversion tool)" ON)
+option(BUILD_PROJ    "Build proj (cartographic projection tool : latlong <-> projected coordinates)" ON)
 
 if(NOT MSVC)
   if (NOT APPLE)
@@ -33,11 +34,6 @@ if(BUILD_CS2CS)
  set(BIN_TARGETS ${BIN_TARGETS} cs2cs)
 endif(BUILD_CS2CS)
 
-if(BUILD_PROJ)
-  include(bin_proj.cmake)
-  set(BIN_TARGETS ${BIN_TARGETS} binproj)
-endif(BUILD_PROJ)
-
 if(BUILD_GEOD)
   include(bin_geod.cmake)
   include(bin_geodtest.cmake)
@@ -49,6 +45,16 @@ if(BUILD_NAD2BIN)
   set(BIN_TARGETS ${BIN_TARGETS} nad2bin)
 endif(BUILD_NAD2BIN)
 
+if(BUILD_PROJ)
+  include(bin_proj.cmake)
+  set(BIN_TARGETS ${BIN_TARGETS} binproj)
+endif(BUILD_PROJ)
+
+if(BUILD_GIE)
+ include(bin_gie.cmake)
+ set(BIN_TARGETS ${BIN_TARGETS} gie)
+endif(BUILD_GIE)
+
 if (MSVC OR CMAKE_CONFIGURATION_TYPES)
   if(BIN_TARGETS)
     # Add _d suffix for your debug versions of the tools
diff --git a/src/Makefile.am b/src/Makefile.am
index 25b7456c..1a1f3270 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -12,7 +12,7 @@ AM_CPPFLAGS =	-DPROJ_LIB=\"$(pkgdatadir)\" \
 include_HEADERS = proj.h proj_api.h projects.h geodesic.h \
 	org_proj4_Projections.h org_proj4_PJ.h
 
-EXTRA_DIST = makefile.vc proj.def bin_cct.cmake bin_cs2cs.cmake \
+EXTRA_DIST = makefile.vc proj.def bin_cct.cmake bin_gie.cmake bin_cs2cs.cmake \
 			 bin_geod.cmake bin_nad2bin.cmake bin_proj.cmake \
 			 lib_proj.cmake CMakeLists.txt bin_geodtest.cmake geodtest.c
 
@@ -21,15 +21,19 @@ cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
 cct_SOURCES = cct.c proj_strtod.c optargpm.h
 nad2bin_SOURCES = nad2bin.c
 geod_SOURCES = geod.c geod_set.c geod_interface.c geod_interface.h
+
+gie_SOURCES = gie.c proj_strtod.c optargpm.h
 multistresstest_SOURCES = multistresstest.c
 test228_SOURCES = test228.c
 geodtest_SOURCES = geodtest.c
 
-proj_LDADD = libproj.la
-cs2cs_LDADD = libproj.la
 cct_LDADD = libproj.la
-nad2bin_LDADD = libproj.la
+cs2cs_LDADD = libproj.la
 geod_LDADD = libproj.la
+proj_LDADD = libproj.la
+nad2bin_LDADD = libproj.la
+
+gie_LDADD = libproj.la
 multistresstest_LDADD = libproj.la @THREAD_LIB@
 test228_LDADD = libproj.la @THREAD_LIB@
 geodtest_LDADD = libproj.la
diff --git a/src/bin_gie.cmake b/src/bin_gie.cmake
new file mode 100644
index 00000000..ca6dde0e
--- /dev/null
+++ b/src/bin_gie.cmake
@@ -0,0 +1,9 @@
+set(GIE_SRC gie.c proj_strtod.c) 
+set(GIE_INCLUDE optargpm.h)
+
+source_group("Source Files\\Bin" FILES ${GIE_SRC})
+
+add_executable(gie ${GIE_SRC} ${GIE_INCLUDE})
+target_link_libraries(gie ${PROJ_LIBRARIES})
+install(TARGETS gie 
+        RUNTIME DESTINATION ${BINDIR})
diff --git a/src/gie.c b/src/gie.c
new file mode 100644
index 00000000..adfc4e0f
--- /dev/null
+++ b/src/gie.c
@@ -0,0 +1,722 @@
+/***********************************************************************
+
+       gie - The Geospatial Integrity Investigation Environment
+
+************************************************************************
+
+The Geospatial Integrity Investigation Environment "gie" is a modest
+regression testing environment for the PROJ.4 transformation library.
+
+Its primary design goal was to be able to replace those thousands of
+lines of regression testing code that are (at time of writing) part
+of PROJ.4, while not requiring any other kind of tooling than the same
+C compiler already employed for compiling the library.
+
+The basic functionality of the gie command language is implemented
+through just 3 command verbs:
+
+OPERATION,     which defines the PROJ.4 operation to test,
+ACCEPT,        which defines the input coordinate to read, and
+EXPECT,        which defines the result to expect.
+
+E.g:
+
+operation  +proj=utm  +zone=32  +ellps=GRS80
+accept     12  55
+expect     691_875.632_14   6_098_907.825_05
+
+Note that gie accepts the underscore ("_") as a thousands separator.
+It is not required (in fact, it is entirely ignored by the input
+routine), but it significantly improves the readability of the very
+long strings of numbers typically required in projected coordinates.
+
+By default, gie considers the EXPECTation met, if the result comes to
+within 0.5 mm of the expected. This default can be changed using the
+TOLERANCE command verb (and yes, I know, linguistically speaking, both
+"operation" and "tolerance" are nouns, not verbs). See the first
+few hundred lines of the "builtins.gie" test file for more details of
+the command verbs available (verbs of both the VERBal and NOUNistic
+persuation).
+
+--
+
+But more importantly than being an acronym for "Geospatial Integrity
+Investigation Environment", gie were also the initials, user id, and
+USGS email address of Gerald Ian Evenden (1935--2016), the geospatial
+visionary, who, already in the 1980s, started what was to become the
+PROJ.4 of today.
+
+Gerald's clear vision was that map projections are *just special
+functions*. Some of them rather complex, most of them of two variables,
+but all of them *just special functions*, and not particularly more
+special than the sin(), cos(), tan(), and hypot() already available in
+the C standard library.
+
+And hence, *they should not be particularly much harder to use*, for a
+programmer, than the sin()s, tan()s and hypot()s so readily available.
+
+Gerald's ingenuity also showed in the implementation of the vision,
+where he devised a highly comprehensible, yet simple, system of key-value
+pairs for parameterising a map projection, and the highly flexible
+PJ struct, storing run-time compiled versions of those key-value pairs,
+hence making a map projection function call, pj_fwd(PJ, point), as easy
+as a traditional function call like hypot(x,y).
+
+While today, we may have more formally well defined metadata systems
+(most prominent the OGC WKT representation), nothing comes close being
+as easily readable ("human compatible") as Gerald's key-value system.
+This system in particular, and the PROJ.4 system in general, was
+Gerald's great gift to anyone using and/or communicating about geodata.
+
+It is only reasonable to name a program keeping an eye on the integrity
+of the PROJ.4 system in honour of Gerald. So in honour, and hopefully
+also in the spirit, of Gerald Ian Evenden (1935--2016), this is the
+Geospatial Integrity Investigation Environment.
+
+************************************************************************
+
+Thomas Knudsen, thokn@sdfe.dk, 2017-10-01/2017-10-08
+
+************************************************************************
+
+* Copyright (c) 2017 Thomas Knudsen
+* Copyright (c) 2017, SDFE
+*
+* Permission is hereby granted, free of charge, to any person obtaining a
+* copy of this software and associated documentation files (the "Software"),
+* to deal in the Software without restriction, including without limitation
+* the rights to use, copy, modify, merge, publish, distribute, sublicense,
+* and/or sell copies of the Software, and to permit persons to whom the
+* Software is furnished to do so, subject to the following conditions:
+*
+* The above copyright notice and this permission notice shall be included
+* in all copies or substantial portions of the Software.
+*
+* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+* DEALINGS IN THE SOFTWARE.
+
+***********************************************************************/
+
+#include "optargpm.h"
+
+#include <proj.h>
+#include "proj_internal.h"
+#include "projects.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+
+#include <string.h>
+#include <ctype.h>
+
+#include <math.h>
+#include <errno.h>
+
+
+
+/* from proj_strtod.c */
+double proj_strtod(const char *str, char **endptr);
+double proj_atof(const char *str);
+
+static char *column (char *buf, int n);
+int   main(int argc, char **argv);
+
+static int   process_file (char *fname);    
+static int   errmsg (int errlev, char *msg, ...);
+static int   get_inp (FILE *f, char *inp, int size);
+static int   get_cmnd (char *inp, char *cmnd, int len);
+static char *get_args (char *inp);
+static int   dispatch (char *cmnd, char *args); 
+
+
+
+#define SKIP -1
+
+typedef struct {
+    char operation[10000];
+    PJ *P;
+    PJ_COORD a, b, c, e;
+    PJ_DIRECTION dir;
+    int verbosity;
+    int nargs;
+    int op_id;
+    int op_ok,    op_ko;
+    int total_ok, total_ko;
+    int grand_ok, grand_ko;
+    double tolerance;
+    char *curr_file;
+    FILE *fout;
+} gie_ctx;
+
+gie_ctx T = {{""}, 0, {{0,0,0,0}}, {{0,0,0,0}}, {{0,0,0,0}}, {{0,0,0,0}}, PJ_FWD, 1,   0,   0,0,0,0,0,0,0,  0.0005, 0, 0};
+
+OPTARGS *o;
+
+
+
+size_t tol_lineno = 0;
+size_t lineno = 0;
+size_t level  = 0;
+char delim[] = {"-------------------------------------------------------------------------------\n"};
+char DELIM[] = {"===============================================================================\n"};
+
+
+#define CMDLEN 25000
+
+int nfiles = 0;
+
+
+static const char usage[] = {
+    "--------------------------------------------------------------------------------\n"
+    "Usage: %s [-options]... infile...\n"
+    "--------------------------------------------------------------------------------\n"
+    "Options:\n"
+    "--------------------------------------------------------------------------------\n"
+    "    -h                Help: print this usage information\n"
+    "    -o /path/to/file  Specify output file name\n"
+    "    -v                Verbose: Provide non-essential informational output.\n"
+    "                      Repeat -v for more verbosity (e.g. -vv)\n"
+    "    -q                Quiet: Opposite of verbose. In quiet mode not even errors\n"
+    "                      are reported. Only interaction is through the return code\n"
+    "                      (0 on success, non-zero indicates number of FAILED tests)\n"
+    "--------------------------------------------------------------------------------\n"
+    "Long Options:\n"
+    "--------------------------------------------------------------------------------\n"
+    "    --output          Alias for -o\n"
+    "    --verbose         Alias for -v\n"
+    "    --help            Alias for -h\n"
+    "--------------------------------------------------------------------------------\n"
+    "Examples:\n"
+    "--------------------------------------------------------------------------------\n"
+    "1. Run all tests in file \"corner-cases.gie\", providing much extra information\n"
+    "       gie -vvvv corner-cases.gie\n"
+    "2. Run all tests in files \"foo\" and \"bar\", providing info on failures only\n"
+    "       gie foo bar\n"
+    "--------------------------------------------------------------------------------\n"    
+};
+
+int main (int argc, char **argv) {
+    int  i;
+    const char *longflags[]  = {"v=verbose", "q=quiet", "h=help", 0};
+    const char *longkeys[]   = {"o=output", 0};
+    
+    o = opt_parse (argc, argv, "hvq", "o", longflags, longkeys);
+    if (0==o)
+        return 0;
+
+    if (opt_given (o, "h")) {
+        printf (usage, o->progname);
+        return 0;
+    }
+
+    
+
+    T.verbosity = opt_given (o, "q");
+    if (T.verbosity)
+        T.verbosity = -1;
+    if (T.verbosity != -1)
+        T.verbosity = opt_given (o, "v") + 1;
+
+    T.fout = stdout;
+    if (opt_given (o, "o"))
+        T.fout = fopen (opt_arg (o, "output"), "rt");
+    if (0==T.fout) {
+        fprintf (stderr, "%s: Cannot open '%s' for output\n", o->progname, opt_arg (o, "output"));
+        free (o);
+        return 1;
+    }
+
+    if (0==o->fargc) {
+        if (T.verbosity==-1)
+            return -1;
+        fprintf (T.fout, "Nothing to do\n");
+        return 0;
+    }
+
+    for (i = 0;  i < o->fargc; i++)
+        process_file (o->fargv[i]);
+
+    if (T.verbosity > 0) {
+        if (o->fargc > 1)
+        fprintf (T.fout, "%sGrand total: %d. Success: %d, Failure: %d\n", delim, T.grand_ok+T.grand_ko, T.grand_ok, T.grand_ko);
+        printf (delim);
+    }
+    else
+        if (T.grand_ko)
+            fprintf (T.fout, "Failures: %d", T.grand_ko);
+
+    if (stdout != T.fout)
+        fclose (T.fout);
+
+    free (o);
+    return T.grand_ko;
+}
+
+
+
+
+    
+static int process_file (char *fname) {
+    FILE *f;
+    char inp[CMDLEN];
+    char cmnd[1000];
+    char *args;
+
+    lineno = level = 0;
+    T.op_ok = T.total_ok = 0;
+    T.op_ko = T.total_ko = 0;
+
+    f = fopen (fname, "rt");
+    if (0==f) {
+        if (T.verbosity > 0) {
+            fprintf (T.fout, "%sCannot open spec'd input file '%s' - bye!\n", delim, fname);
+            return 2;
+        }
+        errmsg (2, "Cannot open spec'd input file '%s' - bye!\n", fname);
+    }
+    if (T.verbosity > 0)
+        fprintf (T.fout, "%sReading file '%s'\n", delim, fname);
+    T.curr_file = fname;
+    while (get_inp(f, inp, CMDLEN)) {
+        int len;
+
+        if (feof(f))
+            break;
+        len = get_cmnd (inp, cmnd, 1000);
+        if (len>=999) {
+            errmsg (2, "Command verb too long: '%s' - bye!\n", cmnd);
+            proj_destroy (T.P);
+            T.P = 0;
+            return 0;
+            
+        }
+        args = get_args (inp);
+        if (SKIP==dispatch (cmnd, args))
+            return proj_destroy (T.P), T.P = 0, 0;
+    }
+    fclose (f);
+
+    T.grand_ok += T.total_ok;
+    T.grand_ko += T.total_ko;
+    if (T.verbosity > 0)
+    fprintf (T.fout, "%stotal: %2d tests succeeded,  %2d tests %s\n", delim, T.total_ok, T.total_ko, T.total_ko? "FAILED!": "failed.");
+    
+    if (level==0)
+        return errmsg (-3, "File '%s':Missing 'BEGIN' cmnd - bye!\n", fname);
+    if (level && level%2)
+        return errmsg (-4, "File '%s':Missing 'END' cmnd - bye!\n", fname);
+    return 0;
+}
+
+
+
+
+
+
+/* return a pointer to the n'th column of buf or a pointer to the terminating 0 if less than n */
+static char *column (char *buf, int n) {
+    int i;
+    if (n <= 0)
+        return buf;
+    for (i = 0;  i < n;  i++) {
+        while (isspace(*buf))
+            buf++;
+        if (i == n - 1)
+            break;
+        while ((0 != *buf) && !isspace(*buf))
+            buf++;
+    }
+    return buf;
+}
+
+
+
+
+static int banner (char *args) {
+    char dots[] = {"..."}, nodots[] = {""}, *thedots = nodots;
+    if (T.total_ko > 0 && T.op_ko==0)
+        printf ("\n\n");
+    if (strlen(args) > 70)
+        thedots = dots;
+        fprintf (T.fout, "%s%-70.70s%s\n", delim, args, thedots);
+    return 0;
+}
+
+
+
+
+
+
+static int tolerance (char *args) {
+    char *endp = args;
+    T.tolerance = proj_strtod (endp, &endp);
+    if (HUGE_VAL==T.tolerance) {
+        T.tolerance = 0.0005;
+        return 1;
+    }
+    while (isspace (*endp))
+        endp++;
+
+    if (0==strcmp(endp, "km"))
+        T.tolerance *= 1000;
+    else if (0==strcmp(endp, "m"))
+        T.tolerance *= 1;
+    else if (0==strcmp(endp, "dm"))
+        T.tolerance /= 10;
+    else if (0==strcmp(endp, "cm"))
+        T.tolerance /= 100;
+    else if (0==strcmp(endp, "mm"))
+        T.tolerance /= 1000;
+    else if (0==strcmp(endp, "um"))
+        T.tolerance /= 1e6;
+    else if (0==strcmp(endp, "nm"))
+        T.tolerance /= 1e9;
+    else
+        T.tolerance /= 1000;  /* If no unit, assume mm */
+    return 0;
+}
+
+
+static int direction (char *args) {
+    char *endp = args;
+    while (isspace (*endp))
+        endp++;
+    switch (*endp) {
+        case 'F':
+        case 'f':
+        T.dir = PJ_FWD;
+            break;
+        case 'I':
+        case 'i':
+        case 'R':
+        case 'r':
+            T.dir = PJ_INV;
+            break;
+        default:
+            return 1;                
+    }
+    return 0;
+}
+
+
+
+
+static void finish_previous_operation () {
+    if (T.verbosity > 1 && T.op_id > 1 && T.op_ok+T.op_ko)
+        fprintf (T.fout, "%s     %d tests succeeded,  %d tests %s\n", delim, T.op_ok, T.op_ko, T.op_ko? "FAILED!": "failed.");
+}
+
+static int operation (char *args) {
+    T.op_id++;
+    strcpy (&(T.operation[0]), args);
+    if (T.verbosity > 1) {
+        finish_previous_operation (args);
+        banner (args);
+    }
+    /* if (0==T.op_ko)
+        printf ("%d\n", (int) tol_lineno); */
+    T.op_ok = 0;
+    T.op_ko = 0;
+    
+    direction ("forward");
+    tolerance ("0.5");
+    
+    if (T.P)
+        proj_destroy (T.P);
+    T.P = proj_create (0, args);
+    if (0==T.P)
+        return errmsg(3, "Invalid operation definition!\n");
+    return 0;
+}
+
+
+static PJ_COORD torad_if_needed (PJ *P, PJ_DIRECTION dir, PJ_COORD a) {
+    enum pj_io_units u = P->left;
+    PJ_COORD c;
+    if (dir==PJ_INV)
+        u = P->right;
+    if (u==PJ_IO_UNITS_CLASSIC || u==PJ_IO_UNITS_METERS)
+        return a;
+    
+    if (u==PJ_IO_UNITS_RADIANS) {
+        c.lpz.lam = proj_torad (T.a.lpz.lam);
+        c.lpz.phi = proj_torad (T.a.lpz.phi);
+    }
+
+    return c;
+}
+
+
+static int accept (char *args) {
+    int n, i;
+    char *endp = args;
+    T.a = proj_coord (0,0,0,0);
+    n = 4;
+    for (i = 0; i < 4; i++) {
+        T.a.v[i] = proj_strtod (endp, &endp);
+        if (HUGE_VAL==T.a.v[i]) {
+            n--;
+            T.a.v[i] = 0;
+        }
+    }
+    T.a = torad_if_needed (T.P, T.dir, T.a);
+    if (T.verbosity > 3)
+        printf ("#  %s", args);
+    return 0;
+}
+
+
+
+static int roundtrip (char *args) {
+    int ntrips;
+    double d, r, ans;
+    char *endp;
+    ans = proj_strtod (args, &endp);
+    ntrips = ans==HUGE_VAL? 100: fabs(ans);
+    d = proj_strtod (endp, &endp);
+    d = d==HUGE_VAL?  T.tolerance:  d / 1000;
+    r = proj_roundtrip (T.P, PJ_FWD, ntrips, T.a);
+    if (r > d) {
+        if (T.verbosity > -1) {
+            if (0==T.op_ko && T.verbosity < 2)
+                banner (T.operation);            
+            fprintf (T.fout, T.op_ko? "     -----\n": delim);
+            fprintf (T.fout, "     FAILURE in %s(%d):\n", opt_strip_path (T.curr_file), (int) lineno);
+            fprintf (T.fout, "     roundtrip deivation: %.3f mm, expected: %.3f mm\n", 1000*r, 1000*d);
+        }
+        T.op_ko++;
+        T.total_ko++;
+    }
+    return 0;
+}
+
+static int expect (char *args) {
+    double d;
+    enum pj_io_units unit;
+    char *endp = args;
+    int i;
+
+    T.e     =  proj_coord (0,0,0,0);
+    T.b     =  proj_coord (0,0,0,0);
+    T.nargs = 4;
+    for (i = 0; i < 4; i++) {
+        T.e.v[i] = proj_strtod (endp, &endp);
+        if (HUGE_VAL==T.e.v[i]) {
+            T.nargs--;
+            T.e.v[i] = 0;
+        }
+    }
+    T.e = torad_if_needed (T.P, T.dir==PJ_FWD? PJ_INV:PJ_FWD, T.e);
+    
+    T.b = proj_trans_coord (T.P, T.dir, T.a);
+    T.b = torad_if_needed (T.P, T.dir==PJ_FWD? PJ_INV:PJ_FWD, T.b);
+    
+    if (T.nargs < 2) {
+        T.op_ko++;
+        T.total_ko++;
+        if (T.verbosity > -1) {
+            if (0==T.op_ko && T.verbosity < 2)
+                banner (T.operation);            
+            fprintf (T.fout, T.op_ko? "     -----\n": delim);
+            fprintf (T.fout, "     FAILURE in %s(%d):\n     Too few args: %s\n", opt_strip_path (T.curr_file), (int) lineno, args);
+        }
+     return 1;
+    }
+
+    unit = T.dir==PJ_FWD? T.P->right: T.P->left;
+    if (PJ_IO_UNITS_CLASSIC==unit)
+        unit = PJ_IO_UNITS_METERS;
+
+    if (unit==PJ_IO_UNITS_RADIANS)
+        d = proj_lp_dist (T.P, T.b.lp, T.e.lp);
+    else
+        d = proj_xyz_dist (T.b.xyz, T.e.xyz);
+
+    if (d > T.tolerance) {
+        if (d > 10)
+            d = 9.999999;
+        if (T.verbosity > -1) {
+            if (0==T.op_ko && T.verbosity < 2)
+                banner (T.operation);
+            fprintf (T.fout, T.op_ko? "     -----\n": delim);
+                
+            fprintf (T.fout, "     FAILURE in %s(%d):\n", opt_strip_path (T.curr_file), (int) lineno);
+            fprintf (T.fout, "     expected: %s\n", args);
+            fprintf (T.fout, "     got:      %.9f   %.9f", T.b.xy.x,  T.b.xy.y);
+            if (T.nargs > 2)
+                fprintf (T.fout, "   %.9f", T.b.xyz.z);
+            if (T.nargs > 3)
+                fprintf (T.fout, "   %.9f", T.b.xyzt.t);
+            fprintf (T.fout, "\n");
+            fprintf (T.fout, "     deviation:  %.3f mm,  expected:  %.3f mm\n", 1000*d, 1000*T.tolerance);
+        }
+        T.op_ko++;
+        T.total_ko++;
+    }
+    else {
+        T.op_ok++;
+        T.total_ok++;
+    }
+    return 0;
+}
+
+
+
+
+
+
+static int verbose (char *args) {
+    int i = proj_atof (args);
+
+    /* if -q/--quiet flag has been given, we do nothing */
+    if (T.verbosity < 0)
+        return 0;
+
+    if (strlen (args))
+        T.verbosity = i;
+    else
+        T.verbosity++;
+    return 0;
+}
+
+static int comment (char *args) {
+    (void) args;
+    return 0;
+}
+
+
+static int echo (char *args) {
+    fprintf (T.fout, "%s\n", args);
+    return 0;
+}
+
+
+
+static int dispatch (char *cmnd, char *args) {
+    if  (0==level%2) {
+        if (0==strcmp (cmnd, "BEGIN"))
+           level++;
+        return 0;
+    }
+    if  (0==stricmp (cmnd, "OPERATION")) return  operation (args);
+    if  (0==stricmp (cmnd, "ACCEPT"))    return  accept    (args);
+    if  (0==stricmp (cmnd, "EXPECT"))    return  expect    (args);
+    if  (0==stricmp (cmnd, "ROUNDTRIP")) return  roundtrip (args);
+    if  (0==stricmp (cmnd, "BANNER"))    return  banner    (args);
+    if  (0==stricmp (cmnd, "VERBOSE"))   return  verbose   (args);
+    if  (0==stricmp (cmnd, "DIRECTION")) return  direction (args);
+    if  (0==stricmp (cmnd, "TOLERANCE")) return  tolerance (args);
+    if  (0==stricmp (cmnd, "ECHO"))      return  echo      (args);
+    if  (0==strcmp  (cmnd, "END"))       return          finish_previous_operation (args), level++, 0;
+    if  ('#'==cmnd[0])                   return  comment   (args);
+    return 0;
+}
+
+
+
+
+
+
+
+
+static int errmsg (int errlev, char *msg, ...) {
+    va_list args;
+    va_start(args, msg);
+    vfprintf(stdout, msg, args);
+    va_end(args);
+    if (errlev)
+        errno = errlev;
+    return errlev;
+}
+
+#define skipspace(f, c)                           \
+    do {                                          \
+        while (isspace (c=fgetc(f)) && !feof(f)){ \
+            if ('\n'==c) lineno++;                \
+        }                                         \
+        if (feof(f))                              \
+            break;                                \
+    } while (ungetc(c, f), 0)
+
+#define skipline(f, c)                            \
+    do {                                          \
+        while ((c=fgetc(f)) && !feof(f)) {        \
+            if ((c=='\r') || (c=='\n'))           \
+                break;                            \
+        }                                         \
+        skipspace (f, c);                         \
+    } while (0)
+
+
+/* skip whitespace at continuation line */
+#define continuation(f, buf, c)                   \
+    if ((c=='\r')||(c=='\n')) {                   \
+        if (c=='\n') lineno++;                    \
+        next--;                                   \
+        while (isspace (c=fgetc(f)) && !feof(f)); \
+    }
+                
+static int get_inp (FILE *f, char *inp, int size) {
+    char *next;
+    int c, esc;
+    char *endp = inp + size - 2;
+
+    skipspace (f, c);
+
+    for (c = esc = 0, next = inp; !feof(f); ) {
+        c = fgetc(f);
+        if (esc) {
+            continuation (f, next, c);
+            esc = 0;
+            /* handle escape sequences here */
+            switch (c) {
+                case '\\': c = '\\'; break;
+                default: (void) c;
+            }
+        }
+        if (c=='\r')
+            break;
+        if (c=='\n') {
+            lineno++;
+            break;
+        }
+            
+        *next++ = c;
+        if ('\\'==c)
+            esc = 1;
+        if (feof(f) || (next==endp))
+            break;
+    }
+    *(next) = 0;
+    return strlen(inp);
+}
+
+static int get_cmnd (char *inp, char *cmnd, int len) {
+    cmnd[0] = 0;
+    while (isspace(*inp++));
+    inp--;
+    while (len-- && !isspace(*inp) && *inp)
+        *cmnd++ = *inp++;
+    *cmnd = 0;
+    return len;
+}
+
+static char *get_args (char *inp) {
+    char *args = inp;
+    while (isspace(*args++))
+        if (0==*args)
+            return args;
+    while (!isspace(*++args))
+        if (0==*args)
+            return args;
+    while (isspace(*args++))
+        if (0==*args)
+            return args;
+    return --args;
+}
diff --git a/src/makefile.vc b/src/makefile.vc
index ac5acf15..fdf03bd3 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -71,12 +71,15 @@ PROJEXE_OBJ	= proj.obj gen_cheb.obj p_series.obj emess.obj
 CS2CSEXE_OBJ	= cs2cs.obj gen_cheb.obj p_series.obj emess.obj
 GEODEXE_OBJ	= geod.obj geod_set.obj geod_interface.obj emess.obj
 CCTEXE_OBJ = cct.obj proj_strtod.obj
+GIEEXE_OBJ = gie.obj proj_strtod.obj
+
 MULTISTRESSTEST_OBJ = multistresstest.obj
 PROJ_DLL 	= proj$(VERSION).dll
 PROJ_EXE    = proj.exe
 CS2CS_EXE   = cs2cs.exe
 GEOD_EXE    = geod.exe
 CCT_EXE     = cct.exe
+GIE_EXE     = gie.exe
 NAD2BIN_EXE = nad2bin.exe
 MULTISTRESSTEST_EXE = multistresstest.exe
 
@@ -85,7 +88,7 @@ CFLAGS	=	/nologo -I. -DPROJ_LIB=\"$(PROJ_LIB_DIR)\" \
 
 default:	all
 
-all: proj.lib $(PROJ_EXE) $(CS2CS_EXE) $(GEOD_EXE) $(CCT_EXE) $(NAD2BIN_EXE)
+all: proj.lib $(PROJ_EXE) $(CS2CS_EXE) $(GEOD_EXE) $(CCT_EXE) $(GIE_EXE) $(NAD2BIN_EXE)
 
 proj.lib:	$(LIBOBJ)
 	if exist proj.lib del proj.lib
@@ -114,6 +117,10 @@ $(CCT_EXE):	$(CCTEXE_OBJ) $(EXE_PROJ)
 	cl $(CCTEXE_OBJ) $(EXE_PROJ)
 	if exist $(CCT_EXE).manifest mt -manifest $(CCT_EXE).manifest -outputresource:$(CCT_EXE);1
 
+$(GIE_EXE):	$(GIEEXE_OBJ) $(EXE_PROJ)
+	cl $(GIEEXE_OBJ) $(EXE_PROJ)
+	if exist $(GIE_EXE).manifest mt -manifest $(GIE_EXE).manifest -outputresource:$(GIE_EXE);1
+
 $(NAD2BIN_EXE):	nad2bin.obj emess.obj $(EXE_PROJ)
 	cl nad2bin.obj emess.obj $(EXE_PROJ)
 
diff --git a/src/optargpm.h b/src/optargpm.h
index 6be2c9ef..4933cd2d 100644
--- a/src/optargpm.h
+++ b/src/optargpm.h
@@ -372,7 +372,8 @@ static int opt_ordinal (OPTARGS *opt, char *option) {
         }
 
     }
-
+    /* kill some potential compiler warnings about unused functions */
+    (void) opt_eof (0);
     return 0;
 }
 
@@ -397,12 +398,23 @@ char *opt_arg (OPTARGS *opt, char *option) {
     return opt->optarg[ordinal];
 }
 
+char *opt_strip_path (char *full_name) {
+    char *last_path_delim, *stripped_name = full_name;
+
+    last_path_delim = strrchr (stripped_name, '\\');
+    if (last_path_delim > stripped_name)
+        stripped_name = last_path_delim + 1;
+
+    last_path_delim = strrchr (stripped_name, '/');
+    if (last_path_delim > stripped_name)
+        stripped_name = last_path_delim + 1;
+    return stripped_name;
+}
 
 /* split command line options into options/flags ("-" style), projdefs ("+" style) and input file args */
 OPTARGS *opt_parse (int argc, char **argv, const char *flags, const char *keys, const char **longflags, const char **longkeys) {
     int i, j;
     OPTARGS *o;
-    char *last_path_delim;
     
     o = (OPTARGS *) calloc (1, sizeof(OPTARGS));
     if (0==o)
@@ -410,14 +422,15 @@ OPTARGS *opt_parse (int argc, char **argv, const char *flags, const char *keys,
 
     o->argc = argc;
     o->argv = argv;
-    o->progname = argv[0];
+    o->progname = opt_strip_path (argv[0]);
+/*    o->progname = argv[0];
     last_path_delim = strrchr (argv[0], '\\');
     if (last_path_delim > o->progname)
         o->progname = last_path_delim;
     last_path_delim = strrchr (argv[0], '/');
     if (last_path_delim > o->progname)
         o->progname = last_path_delim;
-        
+*/        
 
     /* Reset all flags */
     for (i = 0;  i < (int) strlen (flags);  i++)
diff --git a/src/proj_strtod.c b/src/proj_strtod.c
index e9942e5c..c771f2a6 100644
--- a/src/proj_strtod.c
+++ b/src/proj_strtod.c
@@ -86,6 +86,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-01-17/2017-09-18
 ***********************************************************************/
 
 
+#include <string.h> /* for strchr */
 #include <errno.h>
 #include <ctype.h>
 #include <float.h>  /* for HUGE_VAL */
@@ -96,14 +97,16 @@ double proj_atof(const char *str);
     
 
 double proj_strtod(const char *str, char **endptr) {
-    double number = 0;
+    double number = 0, integral_part = 0;
     int exponent = 0;
+    int fraction_is_nonzero = 0;
     int sign = 0;
     char *p = (char *) str;
     int n = 0;
-    int num_digits_total = 0;
-    int num_digits_after_comma = 0;
-
+    int num_digits_total        = 0;
+    int num_digits_after_comma  = 0;
+    int num_prefixed_zeros      = 0;
+    
     if (0==str) {
         errno = EFAULT;
         if (endptr)
@@ -123,10 +126,10 @@ double proj_strtod(const char *str, char **endptr) {
         return HUGE_VAL;
     }
 
-    /* Then handle optional prefixed sign */
+    /* Then handle optional prefixed sign and skip prefix zeros */
     switch (*p) {
         case '-':
-            sign = -1,  p++;  break;
+            sign = -1,  p++;   break;
         case '+':
             sign =  1,  p++;  break;
         default:
@@ -137,7 +140,11 @@ double proj_strtod(const char *str, char **endptr) {
             errno = EINVAL;
             return HUGE_VAL;
     }
-
+    
+    /* skip prefixed zeros */
+    while ('0'==*p || '_'==*p)
+        p++;
+    
     /* Now expect a (potentially zero-length) string of digits */
     while (isdigit(*p) || ('_'==*p)) {
         if ('_'==*p) {
@@ -148,23 +155,45 @@ double proj_strtod(const char *str, char **endptr) {
         p++;
         num_digits_total++;
     }
-        
+    integral_part = number;
+    
     /* Do we have a fractional part? */
     if ('.'==*p) {
         p++;
+        
+        /* keep on skipping prefixed zeros (i.e. allow writing 1e-20 */
+        /* as 0.00000000000000000001 without losing precision) */
+        if (0==integral_part)
+            while ('0'==*p || '_'==*p) {
+                if ('0'==*p)
+                    num_prefixed_zeros++;
+                p++;
+            }
 
+        /* if the next character is nonnumeric, we have reached the end */
+        if (0==strchr ("0123456789eE+-", *p))
+            return integral_part;
+    
         while (isdigit(*p) || '_'==*p) {
-            if ('_'==*p) {
+            /* Don't let pathologically long fractions destroy precision */
+            if ('_'==*p || num_digits_total > 17) {
                 p++;
                 continue;
             }
+
             number = number * 10. + (*p - '0');
+            if (*p!='0')
+                fraction_is_nonzero = 1;
             p++;
             num_digits_total++;
             num_digits_after_comma++;
         }
 
-        exponent = -num_digits_after_comma;
+        /* Avoid having long zero-tails (4321.000...000) destroy precision */
+        if (fraction_is_nonzero)
+            exponent = -(num_digits_after_comma + num_prefixed_zeros);
+        else
+            number = integral_part;
     }
 
     /* non-digit */
@@ -221,8 +250,11 @@ double proj_strtod(const char *str, char **endptr) {
       return HUGE_VAL;
     }
 
-    number *= pow (10, exponent);
+    /* on some platforms pow() is very slow - so don't call it if exponent==0 */
+    if (exponent)
+        number *= pow (10, exponent);
 
+    /* Did we run into an infinity? */
     if (fabs(number) > DBL_MAX)
         errno = ERANGE;