path: root/src/PJ_horner.c

#define PJ_LIB__
#include <proj.h>
#include <projects.h>
#include <assert.h>
#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 */

/***********************************************************************

    Interfacing to a classic piece of geodetic software

************************************************************************

    gen_pol is a highly efficient, classic implementation of a generic
    2D Horner's Scheme polynomial evaluation routine by Knud Poder and
    Karsten Engsager, originating in the vivid geodetic environment at
    what was then (1960-ish) the Danish Geodetic Institute.

    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
    double pointer) is compelling and "true to the code history":

    It has a beautiful classical 1960s ring to it, not unlike the
    original fft implementations, which revolutionized spectral
    analysis in twenty lines of code.

    The Poder coding sound, as classic 1960s as Phil Spector's Wall
    of Sound, is beautiful and inimitable.

	On the other hand: For the uninitiated, the gen_pol code is hard
    to follow, despite being compact.

    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
	more modern (or at least 1990s) look, introducing a "double 2D
	polynomial" data type, HORNER.

    Despite introducing a new data type for handling the polynomial
    coefficients, great care has been taken to keep the coefficient
    array organization identical to that of gen_pol.

    Hence, on one hand, the HORNER data type helps improving the
    long term maintainability of the code by making the data
    organization more mentally accessible.

    On the other hand, it allows us to preserve the business end of
    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
	around 1960, and Karsten Engsager, starting around 1970. It was
    originally written in Algol 60, later (1980s) reimplemented in C.

    The HORNER data type interface, and the organization as a header
    library was implemented by Thomas Knudsen, starting around 2015.

 ***********************************************************************
 *
 * 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
 *
 * 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.
 *
 *****************************************************************************/


 #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 void    horner_free (HORNER *h);

struct horner {
    int    order;    /* maximum degree of polynomium */
    int    coefs;    /* number of coefficients for each polynomium  */
    double range;    /* radius of the region of validity */

    double *fwd_u;   /* coefficients for the forward transformations */
    double *fwd_v;   /* i.e. latitude/longitude to northing/easting  */

    double *inv_u;   /* coefficients for the inverse transformations */
    double *inv_v;   /* i.e. northing/easting to latitude/longitude  */

    UV *fwd_origin;  /* False longitude/latitude */
    UV *inv_origin;  /* False easting/northing   */
};

/* e.g. degree = 2: a + bx + cy + dxx + eyy + fxy, i.e. 6 coefficients */
#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);
}


static HORNER *horner_alloc (size_t order) {
    /* size_t is unsigned, so we need not check for order > 0 */
    int n = horner_number_of_coefficients(order);

    HORNER *h = horner_calloc (1, sizeof (HORNER));

    if (0==h)
        return 0;

    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));

    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)
        return h;

    /* safe, since all pointers are null-initialized (by calloc) */
    horner_free (h);
    return 0;
}




/**********************************************************************/
static UV horner (const HORNER *transformation, int direction, UV position) {
/***********************************************************************

A reimplementation of the classic Engsager/Poder 2D Horner polynomial
evaluation engine "gen_pol".

This version omits the inimitable Poder "dual autocheck"-machinery,
which here is intended to be implemented at a higher level of the
library: We separate the polynomial evaluation from the quality
control (which, given the limited MTBF for "computing machinery",
typical when Knud Poder invented the dual autocheck method,
was not defensible at that time).

Another difference from the original version is that we return the
result on the stack, rather than accepting pointers to result variables
as input. This results in code that is easy to read:

            projected  = horner (s34j,  1, geographic);
            geographic = horner (s34j, -1, projected );

and experiments have shown that on contemporary architectures, the time
taken for returning even comparatively large objects on the stack (and
the UV is not that large - typically only 16 bytes) is negligibly
different from passing two pointers (i.e. typically also 16 bytes) the
other way.

The polynomium has the form:

P = sum (i = [0 : order])
        sum (j = [0 : order - i])
            pow(par_1, i) * pow(par_2, j) * coef(index(order, i, j))

For numerical stability, the summation is carried out backwards,
summing the tiny high order elements first.

***********************************************************************/

    /* These variable names follow the Engsager/Poder  implementation */
    int     sz;              /* Number of coefficients per polynomial */
    double *tcx, *tcy;                        /* Coefficient pointers */
    double  range; /* Equivalent to the gen_pol's FLOATLIMIT constant */
    double  n, e;
    UV uv_error;
    uv_error.u = uv_error.v = HUGE_VAL;

    if (0==transformation)
        return uv_error;

    /* Check for valid value of direction (-1, 0, 1) */
    switch (direction) {
        case 0:    /*  no-op  */
            return position;
        case  1:   /* forward */
        case -1:   /* inverse */
            break;
        default:   /* invalid */
            errno = EINVAL;
            return uv_error;
    }

    /* Prepare for double Horner */
    sz    =  horner_number_of_coefficients(transformation->order);
    range =  transformation->range;


    if (direction==1) {                                   /* forward */
        tcx = transformation->fwd_u + sz;
        tcy = transformation->fwd_v + sz;
        e   = position.u - transformation->fwd_origin->u;
        n   = position.v - transformation->fwd_origin->v;
    } else {                                              /* inverse */
        tcx = transformation->inv_u + sz;
        tcy = transformation->inv_v + 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;
    }

    /* The melody of this block is straight out of the great Engsager/Poder songbook */
    else {
        int g =  transformation->order;
        int r = g, c;
        double u, v, N, E;

        /* Double Horner's scheme: N = n*Cy*e -> yout, E = e*Cx*n -> xout */
        for (N = *--tcy,  E = *--tcx;    r > 0;    r--) {
            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;
        }

        position.u = E;
        position.v = N;
    }

    return position;
}


#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);

}





#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

}


/* 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,

    /*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,

    /*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
  };

#endif /* HORNER_TEST_ORIGINAL_GEN_POL_CODE */


#endif /* ndef HORNER_HEADER_ONLY */

#ifdef __cplusplus
}
#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);
    return point;
}

static PJ_OBS horner_reverse_obs (PJ_OBS point, PJ *P) {
    point.coo.uv = horner ((HORNER *) P->opaque, -1, point.coo.uv);
    return point;
}






static void *horner_freeup (PJ *P) {                        /* Destructor */
    if (0==P)
        return 0;
    if (0==P->opaque)
        return pj_dealloc (P);
    horner_free ((HORNER *) P->opaque);
    return pj_dealloc(P);
}

static void freeup (PJ *P) {
    horner_freeup (P);
    return;
}


static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) {
    char buf[20], *init, *next;
    int i;
    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;

    for (i = 0; i < ncoefs; i++) {
        if (i > 0) {
            if (','!=*next) {
                pj_log_error (P, "Horner: Malformed polynomium set %s. need %d coefs", param, ncoefs);
                return 0;
            }
            init = ++next;
        }
        coefs[i] = pj_strtod (init, &next);
    }
    return 1;
}


/*********************************************************************/
PJ *PROJECTION(horner) {
/*********************************************************************/
    int   degree = 0, n;
    HORNER *Q;
    P->fwdobs  =  horner_forward_obs;
    P->invobs  =  horner_reverse_obs;
    P->fwd3d   =  0;
    P->inv3d   =  0;
    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)");
        return horner_freeup (P);
    }

    Q = horner_alloc (degree);
    P->opaque = (void *) Q;

    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);
    if (0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2))
        return horner_freeup (P);
    if (0==parse_coefs (P, &Q->range, "range", 1))
        Q->range = 500000;

    return P;
}

#ifndef PJ_SELFTEST
/* selftest stub */
int pj_horner_selftest (void) {return 0;}
#else

char tc32_utm32[] = {
    " +proj=horner"
    " +ellps=intl"
    " +range=500000"
    " +fwd_origin=877605.269066,6125810.306769"
    " +inv_origin=877605.760036,6125811.281773"
    " +deg=4"
    " +fwd_v=6.1258112678e+06,9.9999971567e-01,1.5372750011e-10,5.9300860915e-15,2.2609497633e-19,4.3188227445e-05,2.8225130416e-10,7.8740007114e-16,-1.7453997279e-19,1.6877465415e-10,-1.1234649773e-14,-1.7042333358e-18,-7.9303467953e-15,-5.2906832535e-19,3.9984284847e-19"
    " +fwd_u=8.7760574982e+05,9.9999752475e-01,2.8817299305e-10,5.5641310680e-15,-1.5544700949e-18,-4.1357045890e-05,4.2106213519e-11,2.8525551629e-14,-1.9107771273e-18,3.3615590093e-10,2.4380247154e-14,-2.0241230315e-18,1.2429019719e-15,5.3886155968e-19,-1.0167505000e-18"
    " +inv_v=6.1258103208e+06,1.0000002826e+00,-1.5372762184e-10,-5.9304261011e-15,-2.2612705361e-19,-4.3188331419e-05,-2.8225549995e-10,-7.8529116371e-16,1.7476576773e-19,-1.6875687989e-10,1.1236475299e-14,1.7042518057e-18,7.9300735257e-15,5.2881862699e-19,-3.9990736798e-19"
    " +inv_u=8.7760527928e+05,1.0000024735e+00,-2.8817540032e-10,-5.5627059451e-15,1.5543637570e-18,4.1357152105e-05,-4.2114813612e-11,-2.8523713454e-14,1.9109017837e-18,-3.3616407783e-10,-2.4382678126e-14,2.0245020199e-18,-1.2441377565e-15,-5.3885232238e-19,1.0167203661e-18"
};


int pj_horner_selftest (void) {
    PJ *P;
    PJ_OBS a, b;
    double dist;

    /* The TC32 for "System 45 Bornholm" */
    /*    TC32 -> UTM32"    */
    P = pj_create (tc32_utm32);
    if (0==P)
        return 10;

    a = b = pj_obs_null;
    a.coo.uv.v = 6125305.4245;
    a.coo.uv.u =  878354.8539;

    /* Forward projection */
    b = pj_trans (P, PJ_FWD, a);
    b = pj_trans (P, PJ_INV, b);

    /* Check roundtrip precision for 1 iteration each way */
    dist = pj_roundtrip (P, PJ_FWD, 1, a);
    if (dist > 0.01)
        return 1;
    return 0;
}
#endif