PJ_horner: support for complex polynomia (#482)

* PJ_horner: support for complex polynomia

Add Poder/Engsager dual complex Horner and  corresponding test case.
Removed superfluous test code from original Poder/Engsager gen_pol
implementation.

* Trim code in response to a review by @kbevers

* Clean up a few cases of hard coded constants

enum pj_direction symbols replacing hard coded {-1, 0, 1} integer
constants

1 files changed, 182 insertions, 567 deletions

diff --git a/src/PJ_horner.c b/src/PJ_horner.c
index f39d2c2c..08b64ab9 100644
--- a/src/PJ_horner.c
+++ b/src/PJ_horner.c
@@ -5,10 +5,14 @@
 #include <stddef.h>
 #include <math.h>
 #include <errno.h>
+
 PROJ_HEAD(horner,    "Horner polynomial evaluation");
-#define HORNER_SILENCE
 
-/* The next few hundred lines comprises a direct cut-and-paste from the horner.h header library */
+/* make horner.h interface with proj's memory management */
+#define horner_dealloc(x) pj_dealloc(x)
+#define horner_calloc(n,x) pj_calloc(n,x)
+
+/* The next few hundred lines is mostly cut-and-paste from the horner.h header library */
 
 /***********************************************************************
 
@@ -23,7 +27,7 @@ PROJ_HEAD(horner,    "Horner polynomial evaluation");
 
     The original Poder/Engsager gen_pol implementation (where
     the polynomial degree and two sets of polynomial coefficients
-    are packed together in one compound array, handled via a simple
+    are packed together in one compound array, handled via a plain
     double pointer) is compelling and "true to the code history":
 
     It has a beautiful classical 1960s ring to it, not unlike the
@@ -38,7 +42,7 @@ PROJ_HEAD(horner,    "Horner polynomial evaluation");
 
     Also, since adding metadata and improving maintainability
     of the code are among the implied goals of a current SDFE/DTU Space
-	project, the material in this file introduces a version with a more
+	project, the material in this file introduces a version with a
 	more modern (or at least 1990s) look, introducing a "double 2D
 	polynomial" data type, HORNER.
 
@@ -54,13 +58,6 @@ PROJ_HEAD(horner,    "Horner polynomial evaluation");
     the original gen_pol implementation - although not including the
 	famous "Poder dual autocheck" in all its enigmatic elegance.
 
-	The original code has, however, been included in the conditionally
-	compiled TEST-section.
-
-	This is partially for validation of the revised version, partially
-	to enable more generic enjoyment of an interesting piece of
-	ingenious geodetic code - simplistic and enigmatic at the same time.
-
  **********************************************************************
 
 	The material included here was written by Knud Poder, starting
@@ -72,29 +69,7 @@ PROJ_HEAD(horner,    "Horner polynomial evaluation");
 
  ***********************************************************************
  *
- * The gen_pol routine comes with this legal statement (ISC/OpenBSD):
- *
- * Copyright (c) 2011, National Survey and Cadastre, Denmark
- * (Kort- og Matrikelstyrelsen), kms@kms.dk
- *
- * Permission to use, copy, modify, and/or distribute this software for any
- * purpose with or without fee is hereby granted, provided that the above
- * copyright notice and this permission notice appear in all copies.
- *
- * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
- * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
- * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
- * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
- * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
- * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
- * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
- *
- *
- *******************************************************************************
- *
- * The remaining parts are...
- *
- * Copyright (c) 2016, Thomas Knudsen / Karsten Engsager / SDFE http://www.sdfe.dk
+ * Copyright (c) 2016, SDFE http://www.sdfe.dk / Thomas Knudsen / Karsten Engsager
  *
  * Permission is hereby granted, free of charge, to any person obtaining a
  * copy of this software and associated documentation files (the "Software"),
@@ -117,27 +92,10 @@ PROJ_HEAD(horner,    "Horner polynomial evaluation");
  *****************************************************************************/
 
 
- #ifndef HORNER_H
- #define HORNER_H
- #ifdef __cplusplus
- extern "C" {
- #endif
-
-
-
-#if defined(PROJ_H) || defined(PROJECTS_H)
-#define horner_dealloc(x) pj_dealloc(x)
-#define horner_calloc(n,x) pj_calloc(n,x)
-#else
-#define horner_dealloc(x) free(x)
-#define horner_calloc(n,x) calloc(n,x)
-typedef struct {double u,v;} UV;
-#endif
-
 struct horner;
 typedef struct horner HORNER;
-static UV      horner (const HORNER *transformation, int direction, UV position);
-static HORNER *horner_alloc (size_t order);
+static UV      horner (const HORNER *transformation, enum pj_direction, UV position);
+static HORNER *horner_alloc (size_t order, int complex_polynomia);
 static void    horner_free (HORNER *h);
 
 struct horner {
@@ -151,6 +109,9 @@ struct horner {
     double *inv_u;   /* coefficients for the inverse transformations */
     double *inv_v;   /* i.e. northing/easting to latitude/longitude  */
 
+    double *fwd_c;   /* coefficients for the complex forward transformations */
+    double *inv_c;   /* coefficients for the complex inverse transformations */
+
     UV *fwd_origin;  /* False longitude/latitude */
     UV *inv_origin;  /* False easting/northing   */
 };
@@ -159,111 +120,51 @@ struct horner {
 #define horner_number_of_coefficients(order) \
             (((order + 1)*(order + 2)/2))
 
-static int     horner_degree_u (int order, int index);
-static int     horner_degree_v (int order, int index);
-static int     horner_index (int order, int degree_u, int degree_v);
-
-#ifndef HORNER_HEADER_ONLY
-
-/***************************************************************************/
-static int horner_index (int order, int degree_1, int degree_2) {
-/****************************************************************************
-
-    Returns the index of the polynomial coefficient, C, for the element
-
-              C * pow (c_1, degree_2) * pow (c_2, degree_2),
-
-    given that degree_1 > -1, degree_2 > -1,  degree_1 + degree_2 <= order.
-
-    Otherwise returns -1 and sets errno to EDOM.
-
-    The range of the index is [0 : (order + 1) * (order + 2) / 2 - 1].
-
-    A very important thing to note  is that the order of the coordinates
-    c_1 and c_2 depend on the polynomium:
-
-    For the fwd and inv polynomia for the "u" coordinate,
-    u is first (degree_1), v is second (degree_2).
-    For the fwd and inv polynomia for the "v" coordinate,
-    v is first (degree_1), u is second (degree_2).
-
-****************************************************************************/
-    if ( (degree_1 < 0)  ||  (degree_2 < 0)  ||  (degree_1 + degree_2 > order) ) {
-        errno = EDOM;
-        return -1;
-    }
-
-    return (   horner_number_of_coefficients(order) - 1
-             - (order - degree_1)*(order - degree_1 + 1)/2
-             - (order - degree_1 - degree_2));
-}
-
-#define index_u(h, u, v) horner_index (h->order, u, v)
-#define index_v(h, u, v) horner_index (h->order, v, u)
-
-
-static int horner_degree_u (int order, int index) {
-    int n = horner_number_of_coefficients(order);
-    int i, j;
-    if ((order < 0) || (index >= n)) {
-        errno = EDOM;
-        return -1;
-    }
-    for (i = 0; i <= order; i++)
-        for (j = 0; j <= order - i; j++)
-            if (index == horner_index (order, i, j))
-                return i;
-    return -1;
-}
-
-
-static int horner_degree_v (int order, int index) {
-    int n = horner_number_of_coefficients(order);
-    int i, j;
-    if ((order < 0) || (index >= n)) {
-        errno = EDOM;
-        return -1;
-    }
-    for (i = 0; i <= order; i++)
-        for (j = 0; j <= order - i; j++)
-            if (index == horner_index (order, i, j))
-                return j;
-    return -1;
-}
-
-
-
 
 static void horner_free (HORNER *h) {
     horner_dealloc (h->inv_v);
     horner_dealloc (h->inv_u);
     horner_dealloc (h->fwd_v);
     horner_dealloc (h->fwd_u);
+    horner_dealloc (h->fwd_c);
+    horner_dealloc (h->inv_c);
     horner_dealloc (h);
 }
 
 
-static HORNER *horner_alloc (size_t order) {
+static HORNER *horner_alloc (size_t order, int complex_polynomia) {
     /* size_t is unsigned, so we need not check for order > 0 */
     int n = horner_number_of_coefficients(order);
-
+    int polynomia_ok = 0;
     HORNER *h = horner_calloc (1, sizeof (HORNER));
 
     if (0==h)
         return 0;
 
+    if (complex_polynomia)
+        n = 2*order + 2;
     h->order = order;
     h->coefs = n;
 
-    h->fwd_u = horner_calloc (n, sizeof(double));
-    h->fwd_v = horner_calloc (n, sizeof(double));
-    h->inv_u = horner_calloc (n, sizeof(double));
-    h->inv_v = horner_calloc (n, sizeof(double));
+    if (complex_polynomia) {
+        h->fwd_c = horner_calloc (n, sizeof(double));
+        h->inv_c = horner_calloc (n, sizeof(double));
+        if (h->fwd_c && h->inv_c)
+            polynomia_ok = 1;
+    }
+    else {
+        h->fwd_u = horner_calloc (n, sizeof(double));
+        h->fwd_v = horner_calloc (n, sizeof(double));
+        h->inv_u = horner_calloc (n, sizeof(double));
+        h->inv_v = horner_calloc (n, sizeof(double));
+        if (h->fwd_u && h->fwd_v && h->inv_u && h->inv_v)
+            polynomia_ok = 1;
+    }
 
     h->fwd_origin = horner_calloc (1, sizeof(UV));
     h->inv_origin = horner_calloc (1, sizeof(UV));
 
-    if (h->fwd_u && h->fwd_v && h->inv_u && h->inv_v && h->fwd_origin && h->inv_origin)
+    if (polynomia_ok && h->fwd_origin && h->inv_origin)
         return h;
 
     /* safe, since all pointers are null-initialized (by calloc) */
@@ -275,7 +176,7 @@ static HORNER *horner_alloc (size_t order) {
 
 
 /**********************************************************************/
-static UV horner (const HORNER *transformation, int direction, UV position) {
+static UV horner (const HORNER *transformation, enum pj_direction direction, UV position) {
 /***********************************************************************
 
 A reimplementation of the classic Engsager/Poder 2D Horner polynomial
@@ -325,10 +226,10 @@ summing the tiny high order elements first.
 
     /* Check for valid value of direction (-1, 0, 1) */
     switch (direction) {
-        case 0:    /*  no-op  */
+        case PJ_IDENT:    /*  no-op  */
             return position;
-        case  1:   /* forward */
-        case -1:   /* inverse */
+        case PJ_FWD:   /* forward */
+        case PJ_INV:   /* inverse */
             break;
         default:   /* invalid */
             errno = EINVAL;
@@ -340,7 +241,7 @@ summing the tiny high order elements first.
     range =  transformation->range;
 
 
-    if (direction==1) {                                   /* forward */
+    if (direction==PJ_FWD) {                              /* forward */
         tcx = transformation->fwd_u + sz;
         tcy = transformation->fwd_v + sz;
         e   = position.u - transformation->fwd_origin->u;
@@ -381,446 +282,104 @@ summing the tiny high order elements first.
 }
 
 
-#ifdef HORNER_SILENCE
-/**********************************************************************/
-static int horner_silence (int i) {
-/***********************************************************************
-  useless function that silences coompiler warnings about unused stuff
-***********************************************************************/
-    HORNER *Q;
-    UV uv_error;
-    if (i==0)
-        return i;
-    uv_error.u = uv_error.v = HUGE_VAL;
-    horner(0, 1, uv_error);
-    Q = horner_alloc (2);
-    if (Q)
-        horner_free (Q);
-    if (horner_degree_u (2,1))
-        return horner_degree_v (2,1);
-    return i;
-}
-#endif /* def HORNER_SILENCE */
-
-
-
 
-#ifdef HORNER_TEST_ORIGINAL_GEN_POL_CODE
 
-double fwd_u[] = {1,2,3};
-double fwd_v[] = {4,5,6};
-double inv_u[] = {4,6,5};
-double inv_v[] = {1,3,2};
-UV uv_origin_fwd = {0, 0};
-UV uv_origin_inv = {0, 0};
-HORNER uv = {1, 3, 500000.0, fwd_u, fwd_v, inv_u, inv_v, &uv_origin_fwd, &uv_origin_inv};
 
-void show_test (void);
-void tuut_b_test (void);
-
-
-int main (int argc, char **argv) {
-    int i, j, order = 10;
-    UV res = {1,1};
-
-    if (argc==2)
-        order = atoi (argv[1]);
-
-    printf ("Testing %d combinations\n", horner_number_of_coefficients(order));
-
-    for (i = 0; i <= order; i++)
-        for (j = 0; j <= order - i; j++) {
-            int hh = horner_index (order, i, j);
-            int ii = horner_degree_u (order, hh);
-            int jj = horner_degree_v (order, hh);
-            assert (i==ii);
-            assert (j==jj);
-            printf ("%2.2d  %1d%1du  %1d%1dv\n", hh, i, ii, j, jj);
-        }
-
-    tuut_b_test ();
-
-    puts ("Forward...");
-
-    printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
-    res = horner (&uv, 1, res);
-    printf ("res = {%.4f, %.4f}\n", res.u, res.v);
-    assert ( 6==res.u);  /* fwd_u:  a + bu + cv  =  1 + 2*1 + 3*1  =   6 */
-    assert (15==res.v);  /* fwd_v:  a + bu + cv  =  4 + 5*1 + 6*1  =  15 */
-
-    res.u = 2;
-    printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
-    res = horner (&uv, 1, res);
-    printf ("res = {%.4f, %.4f}\n", res.u, res.v);      /* u = 2,  v = 15 */
-    assert (50==res.u);  /* fwd_u:  a + bu + cv  =  1 + 2*(u=2)  + 3*(v=15)  =  50 */
-    assert (91==res.v);  /* fwd_v:  a + bu + cv  =  4 + 5*(v=15) + 6*(u=2)   =  91 */
-
-    puts ("Inverse...");
-
-    res.u = 1;
-    res.v = 1;
-    printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
-    res = horner (&uv, -1, res);
-    printf ("res = {%.4f, %.4f}\n", res.u, res.v);      /* u = 1,  v = 1 */
-    assert (15==res.u);  /* inv_u:  a + bu + cv  =  4 + 5*1 + 6*1  =  15 */
-    assert ( 6==res.v);  /* inv_v:  a + bu + cv  =  1 + 2*1 + 3*1  =   6 */
-
-    res.u = 2;
-    printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
-    res = horner (&uv, -1, res);
-    printf ("res = {%.4f, %.4f}\n", res.u, res.v);      /* u = 2,  v = 6 */
-    assert (46==res.u);  /* inv_u:  a + bu + cv  =  4 + 6*(u=2) + 5*(v=6)  =  46 */
-    assert (23==res.v);  /* inv_v:  a + bu + cv  =  1 + 2*(u=2) + 3*(v=6)  =  23 */
-
-/*  show_test ();*/
-    return 0;
-}
-
-
-
-
-
-
-/* ttu_n and ttu_e are based on static double C_ttu_b[] */
-static  double ttu_n[] = {
-    /* tc32_ed50 -> utm32_ed50 : Bornholm */
-    /* m_lim_gen:  0.086  red  = 0      OBS   =    852   */
-    /* m = 1.38 cm my_loss =    +2  y_enp = +10.5 */
-    /* m = 1.44 cm mx_loss =    +2  x_enp = +10.4 */
-
-    /*deg  4.0,*/
-    /*Poly NORTH :: e-degree =  0  : n-degree =  4 */
-    /*  0*/   6.1258112678e+06,   9.9999971567e-01,   1.5372750011e-10,
-    /*  3*/   5.9300860915e-15,   2.2609497633e-19,
-
-    /*Poly NORTH :: e-degree =  1  : n-degree =  3 */
-    /*  5*/   4.3188227445e-05,   2.8225130416e-10,   7.8740007114e-16,
-    /*  8*/  -1.7453997279e-19,
-
-    /*Poly NORTH :: e-degree =  2  : n-degree =  2 */
-    /*  9*/   1.6877465415e-10,  -1.1234649773e-14,  -1.7042333358e-18,
-    /*Poly NORTH :: e-degree =  3  : n-degree =  1 */
-    /* 12*/  -7.9303467953e-15,  -5.2906832535e-19,
-    /*Poly NORTH :: e-degree =  4  : n-degree =  0 */
-    /* 14*/   3.9984284847e-19,
-
-    /*tcy    6125810.306769, */
-};
-
-static  double ttu_e[] = {
-    /*Poly EAST  :: n-degree =  0  : e-degree =  4 */
-    /*  0*/   8.7760574982e+05,   9.9999752475e-01,   2.8817299305e-10,
-    /*  3*/   5.5641310680e-15,  -1.5544700949e-18,
-
-    /*Poly EAST  :: n-degree =  1  : e-degree =  3 */
-    /*  5*/  -4.1357045890e-05,   4.2106213519e-11,   2.8525551629e-14,
-    /*  8*/  -1.9107771273e-18,
-
-    /*Poly EAST  :: n-degree =  2  : e-degree =  2 */
-    /*  9*/   3.3615590093e-10,   2.4380247154e-14,  -2.0241230315e-18,
-    /*Poly EAST  :: n-degree =  3  : e-degree =  1 */
-    /* 12*/   1.2429019719e-15,   5.3886155968e-19,
-    /*Poly EAST  :: n-degree =  4  : e-degree =  0 */
-    /* 14*/  -1.0167505000e-18,
-
-    /* tcx     877605.269066 */
-  };
-
-
-/* utt_n and utt_e are based on static  double C_utt_b[] */
-static  double utt_n[] = {
-    /* utm32_ed50 -> tc32_ed50 : Bornholm */
-    /* m_lim_gen:  0.086  red  = 0      OBS   =    852   */
-    /* m = 1.38 cm my_loss =    +2  y_enp = +10.8 */
-    /* m = 1.44 cm mx_loss =    +2  x_enp = +10.7 */
-
-    /*deg  4.0,*/
-    /*Poly NORTH :: e-degree =  0  : n-degree =  4 */
-    /*  0*/   6.1258103208e+06,   1.0000002826e+00,  -1.5372762184e-10,
-    /*  3*/  -5.9304261011e-15,  -2.2612705361e-19,
-
-    /*Poly NORTH :: e-degree =  1  : n-degree =  3 */
-    /*  5*/  -4.3188331419e-05,  -2.8225549995e-10,  -7.8529116371e-16,
-    /*  8*/   1.7476576773e-19,
-
-    /*Poly NORTH :: e-degree =  2  : n-degree =  2 */
-    /*  9*/  -1.6875687989e-10,   1.1236475299e-14,   1.7042518057e-18,
-    /*Poly NORTH :: e-degree =  3  : n-degree =  1 */
-    /* 12*/   7.9300735257e-15,   5.2881862699e-19,
-    /*Poly NORTH :: e-degree =  4  : n-degree =  0 */
-    /* 14*/  -3.9990736798e-19,
-
-    /*tcy    6125811.281773,*/
-};
-
-static  double utt_e[] = {
-    /*Poly EAST  :: n-degree =  0  : e-degree =  4 */
-    /*  0*/   8.7760527928e+05,   1.0000024735e+00,  -2.8817540032e-10,
-    /*  3*/  -5.5627059451e-15,   1.5543637570e-18,
-
-    /*Poly EAST  :: n-degree =  1  : e-degree =  3 */
-    /*  5*/   4.1357152105e-05,  -4.2114813612e-11,  -2.8523713454e-14,
-    /*  8*/   1.9109017837e-18,
-
-    /*Poly EAST  :: n-degree =  2  : e-degree =  2 */
-    /*  9*/  -3.3616407783e-10,  -2.4382678126e-14,   2.0245020199e-18,
-    /*Poly EAST  :: n-degree =  3  : e-degree =  1 */
-    /* 12*/  -1.2441377565e-15,  -5.3885232238e-19,
-    /*Poly EAST  :: n-degree =  4  : e-degree =  0 */
-    /* 14*/   1.0167203661e-18,
-
-    /*tcx     877605.760036 */
-};
-
-UV tuut_b_origin_fwd = {877605.269066, 6125810.306769};
-UV tuut_b_origin_inv = {877605.760036, 6125811.281773};
-HORNER tuut_b = {4, 15, 500000.0,   ttu_e, ttu_n, utt_e, utt_n,   &tuut_b_origin_fwd, &tuut_b_origin_inv};
-
-
-/* Prototype and forward declarations of the material needed for cross-check with original implementation */
-int gen_pol(double *C_f, double *C_r, double N_in, double E_in, double *Nout, double *Eout);
-static  double      C_ttu_b[];
-static  double      C_utt_b[];
-
-
-void gen_pol_roundtrip (double *C_ttu_b, double *C_utt_b, UV fwd) {
-    UV res, hrn;
-    int ret;
-
-    ret = gen_pol(C_ttu_b, C_utt_b, fwd.v, fwd.u, &res.v, &res.u);
-    printf ("\n------\n\n");
-	if (0!=ret) printf ("ret: %d\n", ret);
-    printf ("inp: %11.3f  %11.3f\n", fwd.u, fwd.v);
-    printf ("res: %11.3f  %11.3f\n", res.u, res.v);
-    hrn = horner (&tuut_b, 1, fwd);
-    printf ("hrn: %11.3f  %11.3f\n", hrn.u, hrn.v);
-    assert (hrn.u==res.u);
-    assert (hrn.v==res.v);
-
-    ret = gen_pol(C_utt_b, C_ttu_b, res.v, res.u, &res.v, &res.u);
-    hrn = horner (&tuut_b, -1, hrn);
-    printf ("hrn: %11.3f  %11.3f\n", hrn.u, hrn.v);
-	if (0!=ret) printf ("ret: %d\n", ret);
-    printf ("res: %11.3f  %11.3f\n",   res.u, res.v);
-    assert (hrn.u==res.u);
-    assert (hrn.v==res.v);
-    printf ("inp: %11.3f  %11.3f  (%.3g mm)\n", fwd.u, fwd.v, 1e3*hypot(fwd.u-res.u, fwd.v-res.v));
-}
 
-void tuut_b_test (void) {
-    UV fwd = tuut_b_origin_fwd;
-    UV res;
-
-    puts ("Bornholm");
-    printf ("fwd: %11.3f  %11.3f\n", fwd.u, fwd.v);
-    res = horner (&tuut_b, 1, fwd);
-    printf ("res: %11.3f  %11.3f\n", res.u, res.v);
-    res = horner (&tuut_b, -1, res);
-    printf ("res: %11.3f  %11.3f\n", res.u, res.v);
-    printf ("fwd: %11.3f  %11.3f\n", fwd.u, fwd.v);
-
-    gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
-    fwd.u = 877000;
-    fwd.v = 6125000;
-    gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
-    fwd.u = 800000;
-    fwd.v = 6100000;
-    gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
-    fwd.u = 850000;
-    fwd.v = 6200000;
-    gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
 
+static PJ_OBS horner_forward_obs (PJ_OBS point, PJ *P) {
+    point.coo.uv = horner ((HORNER *) P->opaque, 1, point.coo.uv);
+    return point;
 }
 
-
-
-
-
-#define            FLOATLIMIT      5.0e5
-#define TRF_AREA_ EDOM
-
-/********************************************************************************************/
-int gen_pol(double *C_f, double *C_r, double N_in, double E_in, double *Nout, double *Eout) {
-/*********************************************************************************************
-
-    This is the original Poder/Engsager implementation of gen_pol.
-	It is included here for test-comparison with the horner() routine.
-
-*********************************************************************************************/
-  double             N, E, n, e;
-  double            *Cp, *tcy, *tcx;
-  double             tol = 1.0e-4;
-  int                i;
-  int                g;
-  int                sz;
-  int                r, c;
-  int                res = 0;
-
-  /* Preserve input for reverse check */
-  N  = N_in;
-  E  = E_in;
-  Cp = C_f;
-
-  /* Transformation loop */
-  for (i = -1; i <= 1 && res == 0; ++i)
-    if (i) {
-
-      /* Prepare for double Horner */
-      g   = (int) *Cp;
-      sz  = (g + 1)*(g + 2)/2 + 1;
-      tcy = Cp + sz;
-      tcx = tcy  + sz;
-      /* Double Horner's scheme */
-      /* N = n*Cy*e -> yout, E = e*Cx*n -> xout */
-      n   = N - *tcy;
-      e   = E - *tcx;
-      if ((fabs(n) < FLOATLIMIT) && (fabs(e) < FLOATLIMIT)) {
-
-        for (  r = g, N = *--tcy, E = *--tcx; r >  0; r--) {
-          double u, v;
-          for (c = g, u = *--tcy, v = *--tcx; c >= r; c--) {
-            u = n*u + *--tcy;
-            v = e*v + *--tcx;
-          }
-          N = e*N + u;
-          E = n*E + v;
-        }
-      } else res = TRF_AREA_;
-    }
-    else { /* collect output coord,switch to reverse checking */
-      *Nout = N;
-      *Eout = E;
-      Cp    = C_r;
-    }
-
-  /* tol-check of results*/
-  if (res == 0 && (fabs(N - N_in) < tol && fabs(E - E_in) < tol))
-      return (0);
-  else if (res == 0) res = TRF_AREA_;
-  return(res);
-
-#undef   FLOATLIMIT
-
+static PJ_OBS horner_reverse_obs (PJ_OBS point, PJ *P) {
+    point.coo.uv = horner ((HORNER *) P->opaque, -1, point.coo.uv);
+    return point;
 }
 
 
-/* s45b polynomia */
-  static  double      C_ttu_b[] = {
-    /* tc32_ed50 -> utm32_ed50 : Bornholm */
-    /* m_lim_gen:  0.086  red  = 0      OBS   =    852   */
-    /* m = 1.38 cm my_loss =    +2  y_enp = +10.5 */
-    /* m = 1.44 cm mx_loss =    +2  x_enp = +10.4 */
-
-    /*deg*/  4.0,
-    /*Poly NORTH :: e-degree =  0  : n-degree =  4 */
-    /*  0*/   6.1258112678e+06,   9.9999971567e-01,   1.5372750011e-10,
-    /*  3*/   5.9300860915e-15,   2.2609497633e-19,
-
-    /*Poly NORTH :: e-degree =  1  : n-degree =  3 */
-    /*  5*/   4.3188227445e-05,   2.8225130416e-10,   7.8740007114e-16,
-    /*  8*/  -1.7453997279e-19,
-
-    /*Poly NORTH :: e-degree =  2  : n-degree =  2 */
-    /*  9*/   1.6877465415e-10,  -1.1234649773e-14,  -1.7042333358e-18,
-    /*Poly NORTH :: e-degree =  3  : n-degree =  1 */
-    /* 12*/  -7.9303467953e-15,  -5.2906832535e-19,
-    /*Poly NORTH :: e-degree =  4  : n-degree =  0 */
-    /* 14*/   3.9984284847e-19,
-
-    /*tcy*/    6125810.306769,
-
-    /*Poly EAST  :: n-degree =  0  : e-degree =  4 */
-    /*  0*/   8.7760574982e+05,   9.9999752475e-01,   2.8817299305e-10,
-    /*  3*/   5.5641310680e-15,  -1.5544700949e-18,
-
-    /*Poly EAST  :: n-degree =  1  : e-degree =  3 */
-    /*  5*/  -4.1357045890e-05,   4.2106213519e-11,   2.8525551629e-14,
-    /*  8*/  -1.9107771273e-18,
-
-    /*Poly EAST  :: n-degree =  2  : e-degree =  2 */
-    /*  9*/   3.3615590093e-10,   2.4380247154e-14,  -2.0241230315e-18,
-    /*Poly EAST  :: n-degree =  3  : e-degree =  1 */
-    /* 12*/   1.2429019719e-15,   5.3886155968e-19,
-    /*Poly EAST  :: n-degree =  4  : e-degree =  0 */
-    /* 14*/  -1.0167505000e-18,
-
-    /*tcx*/     877605.269066
-  };
 
-  static  double      C_utt_b[] = {
-    /* utm32_ed50 -> tc32_ed50 : Bornholm */
-    /* m_lim_gen:  0.086  red  = 0      OBS   =    852   */
-    /* m = 1.38 cm my_loss =    +2  y_enp = +10.8 */
-    /* m = 1.44 cm mx_loss =    +2  x_enp = +10.7 */
 
-    /*deg*/  4.0,
-    /*Poly NORTH :: e-degree =  0  : n-degree =  4 */
-    /*  0*/   6.1258103208e+06,   1.0000002826e+00,  -1.5372762184e-10,
-    /*  3*/  -5.9304261011e-15,  -2.2612705361e-19,
-
-    /*Poly NORTH :: e-degree =  1  : n-degree =  3 */
-    /*  5*/  -4.3188331419e-05,  -2.8225549995e-10,  -7.8529116371e-16,
-    /*  8*/   1.7476576773e-19,
+/**********************************************************************/
+static UV complex_horner (const HORNER *transformation, enum pj_direction direction, UV position) {
+/***********************************************************************
 
-    /*Poly NORTH :: e-degree =  2  : n-degree =  2 */
-    /*  9*/  -1.6875687989e-10,   1.1236475299e-14,   1.7042518057e-18,
-    /*Poly NORTH :: e-degree =  3  : n-degree =  1 */
-    /* 12*/   7.9300735257e-15,   5.2881862699e-19,
-    /*Poly NORTH :: e-degree =  4  : n-degree =  0 */
-    /* 14*/  -3.9990736798e-19,
+A reimplementation of a classic Engsager/Poder Horner complex
+polynomial evaluation engine.
 
-    /*tcy*/    6125811.281773,
+***********************************************************************/
 
-    /*Poly EAST  :: n-degree =  0  : e-degree =  4 */
-    /*  0*/   8.7760527928e+05,   1.0000024735e+00,  -2.8817540032e-10,
-    /*  3*/  -5.5627059451e-15,   1.5543637570e-18,
+    /* These variable names follow the Engsager/Poder  implementation */
+    int     sz;                             /* Number of coefficients */
+    double *c, *cb;                           /* Coefficient pointers */
+    double  range; /* Equivalent to the gen_pol's FLOATLIMIT constant */
+    double  n, e, w, N, E;
+    UV uv_error;
+    uv_error.u = uv_error.v = HUGE_VAL;
 
-    /*Poly EAST  :: n-degree =  1  : e-degree =  3 */
-    /*  5*/   4.1357152105e-05,  -4.2114813612e-11,  -2.8523713454e-14,
-    /*  8*/   1.9109017837e-18,
+    if (0==transformation)
+        return uv_error;
 
-    /*Poly EAST  :: n-degree =  2  : e-degree =  2 */
-    /*  9*/  -3.3616407783e-10,  -2.4382678126e-14,   2.0245020199e-18,
-    /*Poly EAST  :: n-degree =  3  : e-degree =  1 */
-    /* 12*/  -1.2441377565e-15,  -5.3885232238e-19,
-    /*Poly EAST  :: n-degree =  4  : e-degree =  0 */
-    /* 14*/   1.0167203661e-18,
+    /* Check for valid value of direction (-1, 0, 1) */
+    switch (direction) {
+        case PJ_IDENT:    /*  no-op  */
+            return position;
+        case PJ_FWD:   /* forward */
+        case PJ_INV:   /* inverse */
+            break;
+        default:   /* invalid */
+            errno = EINVAL;
+            return uv_error;
+    }
 
-    /*tcx*/     877605.760036
-  };
+    /* Prepare for double Horner */
+    sz    =  2*transformation->order + 2;
+    range =  transformation->range;
 
-#endif /* HORNER_TEST_ORIGINAL_GEN_POL_CODE */
+    if (direction==PJ_FWD) {                              /* forward */
+        cb =  transformation->fwd_c;
+        c  =  cb + sz;
+        e  =  position.u - transformation->fwd_origin->u;
+        n  =  position.v - transformation->fwd_origin->v;
+    } else {                                              /* inverse */
+        cb =  transformation->inv_c;
+        c  =  cb + sz;
+        e  =  position.u - transformation->inv_origin->u;
+        n  =  position.v - transformation->inv_origin->v;
+    }
 
+    if ((fabs(n) > range) || (fabs(e) > range)) {
+        errno = EDOM;
+        return uv_error;
+    }
 
-#endif /* ndef HORNER_HEADER_ONLY */
+    /* Everything's set up properly - now do the actual polynomium evaluation */
+    E = *--c;
+    N = *--c;
+    while (c > cb) {
+        w = n*E + e*N + *--c;
+        N = n*N - e*E + *--c;
+        E = w;
+    }
 
-#ifdef __cplusplus
+    position.u = E;
+    position.v = N;
+    return position;
 }
-#endif
 
-#endif /* ndef HORNER_H */
 
 
-
-
-
-
-
-static PJ_OBS horner_forward_obs (PJ_OBS point, PJ *P) {
-    point.coo.uv = horner ((HORNER *) P->opaque, 1, point.coo.uv);
+static PJ_OBS complex_horner_forward_obs (PJ_OBS point, PJ *P) {
+    point.coo.uv = complex_horner ((HORNER *) P->opaque, PJ_FWD, point.coo.uv);
     return point;
 }
 
-static PJ_OBS horner_reverse_obs (PJ_OBS point, PJ *P) {
-    point.coo.uv = horner ((HORNER *) P->opaque, -1, point.coo.uv);
+static PJ_OBS complex_horner_reverse_obs (PJ_OBS point, PJ *P) {
+    point.coo.uv = complex_horner ((HORNER *) P->opaque, PJ_INV, point.coo.uv);
     return point;
 }
 
 
-
-
-
-
 static void *horner_freeup (PJ *P) {                        /* Destructor */
     if (0==P)
         return 0;
@@ -837,13 +396,21 @@ static void freeup (PJ *P) {
 
 
 static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) {
-    char buf[20], *init, *next;
+    char *buf, *init, *next;
     int i;
+
+    buf = pj_calloc (strlen (param) + 2, sizeof(char));
+    if (0==buf) {
+        pj_log_error (P, "Horner: Out of core");
+        return 0;
+    }
+
     sprintf (buf, "t%s", param);
     if (0==pj_param (P->ctx, P->params, buf).i)
         return 0;
     sprintf (buf, "s%s", param);
     init = pj_param(P->ctx, P->params, buf).s;
+    pj_dealloc (buf);
 
     for (i = 0; i < ncoefs; i++) {
         if (i > 0) {
@@ -862,7 +429,7 @@ static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) {
 /*********************************************************************/
 PJ *PROJECTION(horner) {
 /*********************************************************************/
-    int   degree = 0, n;
+    int   degree = 0, n, complex_horner = 0;
     HORNER *Q;
     P->fwdobs  =  horner_forward_obs;
     P->invobs  =  horner_reverse_obs;
@@ -871,30 +438,42 @@ PJ *PROJECTION(horner) {
     P->fwd     =  0;
     P->inv     =  0;
     P->left    =  P->right  =  PJ_IO_UNITS_METERS;
-    /* silence a few compiler warnings */
-    horner_silence (0);
 
     /* Polynomial degree specified? */
     if (pj_param (P->ctx, P->params, "tdeg").i) /* degree specified? */
 		degree = pj_param(P->ctx, P->params, "ideg").i;
     else {
-        pj_log_debug (P, "Horner: Need to specify polynomial degree, (+deg=n)");
+        pj_log_debug (P, "Horner: Must specify polynomial degree, (+deg=n)");
         return horner_freeup (P);
     }
 
-    Q = horner_alloc (degree);
+    if (pj_param (P->ctx, P->params, "tfwd_c").i || pj_param (P->ctx, P->params, "tinv_c").i) /* complex polynomium? */
+		complex_horner = 1;
+
+    Q = horner_alloc (degree, complex_horner);
     P->opaque = (void *) Q;
 
-    n = horner_number_of_coefficients (degree);
+    if (complex_horner) {
+        n = 2*degree + 2;
+        if (0==parse_coefs (P, Q->fwd_c, "fwd_c", n))
+            return horner_freeup (P);
+        if (0==parse_coefs (P, Q->inv_c, "inv_c", n))
+            return horner_freeup (P);
+        P->fwdobs  =  complex_horner_forward_obs;
+        P->invobs  =  complex_horner_reverse_obs;
 
-    if (0==parse_coefs (P, Q->fwd_u, "fwd_u", n))
-        return horner_freeup (P);
-    if (0==parse_coefs (P, Q->fwd_v, "fwd_v", n))
-        return horner_freeup (P);
-    if (0==parse_coefs (P, Q->inv_u, "inv_u", n))
-        return horner_freeup (P);
-    if (0==parse_coefs (P, Q->inv_v, "inv_v", n))
-        return horner_freeup (P);
+    }
+    else {
+        n = horner_number_of_coefficients (degree);
+        if (0==parse_coefs (P, Q->fwd_u, "fwd_u", n))
+            return horner_freeup (P);
+        if (0==parse_coefs (P, Q->fwd_v, "fwd_v", n))
+            return horner_freeup (P);
+        if (0==parse_coefs (P, Q->inv_u, "inv_u", n))
+            return horner_freeup (P);
+        if (0==parse_coefs (P, Q->inv_v, "inv_v", n))
+            return horner_freeup (P);
+    }
 
     if (0==parse_coefs (P, (double *)(Q->fwd_origin), "fwd_origin", 2))
         return horner_freeup (P);
@@ -910,7 +489,6 @@ PJ *PROJECTION(horner) {
 /* selftest stub */
 int pj_horner_selftest (void) {return 0;}
 #else
-
 char tc32_utm32[] = {
     " +proj=horner"
     " +ellps=intl"
@@ -925,13 +503,24 @@ char tc32_utm32[] = {
 };
 
 
+char sb_utm32[] = {
+    " +proj=horner"
+    " +ellps=intl"
+    " +range=500000"
+    " +tolerance=0.0005"
+    " +fwd_origin=4.94690026817276e+05,6.13342113183056e+06"
+    " +inv_origin=6.19480258923588e+05,6.13258568148837e+06"
+    " +deg=3"
+    " +fwd_c=6.13258562111350e+06,6.19480105709997e+05,9.99378966275206e-01,-2.82153291753490e-02,-2.27089979140026e-10,-1.77019590701470e-09,1.08522286274070e-14,2.11430298751604e-15"
+    " +inv_c=6.13342118787027e+06,4.94690181709311e+05,9.99824464710368e-01,2.82279070814774e-02,7.66123542220864e-11,1.78425334628927e-09,-1.05584823306400e-14,-3.32554258683744e-15"
+};
+
 int pj_horner_selftest (void) {
     PJ *P;
-    PJ_OBS a, b;
+    PJ_OBS a, b, c;
     double dist;
 
-    /* The TC32 for "System 45 Bornholm" */
-    /*    TC32 -> UTM32"    */
+    /* Real polynonia relating the technical coordinate system TC32 to "System 45 Bornholm" */
     P = pj_create (tc32_utm32);
     if (0==P)
         return 10;
@@ -940,14 +529,40 @@ int pj_horner_selftest (void) {
     a.coo.uv.v = 6125305.4245;
     a.coo.uv.u =  878354.8539;
 
+    /* Check roundtrip precision for 1 iteration each way, starting in forward direction */
+    dist = pj_roundtrip (P, PJ_FWD, 1, a);
+    if (dist > 0.01)
+        return 1;
+
+    /* The complex polynomial transformation between the "System Storebaelt" and utm32/ed50 */
+    P = pj_create (sb_utm32);
+    if (0==P)
+        return 11;
+
+    /* Test value: utm32_ed50(620000, 6130000) = sb_ed50(495136.8544, 6130821.2945) */
+    a = b = c = pj_obs_null;
+    a.coo.uv.v = 6130821.2945;
+    a.coo.uv.u =  495136.8544;
+    c.coo.uv.v = 6130000.0000;
+    c.coo.uv.u =  620000.0000;
+
     /* Forward projection */
     b = pj_trans (P, PJ_FWD, a);
-    b = pj_trans (P, PJ_INV, b);
+    dist = pj_xy_dist (b.coo.xy, c.coo.xy);
+    if (dist > 0.001)
+        return 2;
+
+    /* Inverse projection */
+    b = pj_trans (P, PJ_INV, c);
+    dist = pj_xy_dist (b.coo.xy, a.coo.xy);
+    if (dist > 0.001)
+        return 3;
 
     /* Check roundtrip precision for 1 iteration each way */
     dist = pj_roundtrip (P, PJ_FWD, 1, a);
     if (dist > 0.01)
-        return 1;
+        return 4;
+
     return 0;
 }
 #endif