EVOLUTION-MANAGER

Edit File: horner.cpp

/*********************************************************************** 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 plain 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 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 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. *********************************************************************** * * 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"), * 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. * *****************************************************************************/ #define PJ_LIB__ #include <errno.h> #include <math.h> #include <stddef.h> #include <stdio.h> #include <string.h> #include "proj.h" #include "proj_internal.h" PROJ_HEAD(horner, "Horner polynomial evaluation"); /* 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) namespace { // anonymous namespace struct horner { int uneg; /* u axis negated? */ int vneg; /* v axis negated? */ 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 */ double *fwd_c; /* coefficients for the complex forward transformations */ double *inv_c; /* coefficients for the complex inverse transformations */ PJ_UV *fwd_origin; /* False longitude/latitude */ PJ_UV *inv_origin; /* False easting/northing */ }; } // anonymous namespace typedef struct horner HORNER; static PJ_UV horner_func (const HORNER *transformation, PJ_DIRECTION direction, PJ_UV position); static HORNER *horner_alloc (size_t order, int complex_polynomia); static void horner_free (HORNER *h); /* 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 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->fwd_origin); horner_dealloc (h->inv_origin); horner_dealloc (h); } static HORNER *horner_alloc (size_t order, int complex_polynomia) { /* size_t is unsigned, so we need not check for order > 0 */ int n = (int)horner_number_of_coefficients(order); int polynomia_ok = 0; HORNER *h = static_cast<HORNER*>(horner_calloc (1, sizeof (HORNER))); if (nullptr==h) return nullptr; if (complex_polynomia) n = 2*(int)order + 2; h->order = (int)order; h->coefs = n; if (complex_polynomia) { h->fwd_c = static_cast<double*>(horner_calloc (n, sizeof(double))); h->inv_c = static_cast<double*>(horner_calloc (n, sizeof(double))); if (h->fwd_c && h->inv_c) polynomia_ok = 1; } else { h->fwd_u = static_cast<double*>(horner_calloc (n, sizeof(double))); h->fwd_v = static_cast<double*>(horner_calloc (n, sizeof(double))); h->inv_u = static_cast<double*>(horner_calloc (n, sizeof(double))); h->inv_v = static_cast<double*>(horner_calloc (n, sizeof(double))); if (h->fwd_u && h->fwd_v && h->inv_u && h->inv_v) polynomia_ok = 1; } h->fwd_origin = static_cast<PJ_UV*>(horner_calloc (1, sizeof(PJ_UV))); h->inv_origin = static_cast<PJ_UV*>(horner_calloc (1, sizeof(PJ_UV))); if (polynomia_ok && h->fwd_origin && h->inv_origin) return h; /* safe, since all pointers are null-initialized (by calloc) */ horner_free (h); return nullptr; } /**********************************************************************/ static PJ_UV horner_func (const HORNER *transformation, PJ_DIRECTION direction, PJ_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; PJ_UV uv_error; uv_error.u = uv_error.v = HUGE_VAL; if (nullptr==transformation) return uv_error; /* 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; } /* Prepare for double Horner */ sz = horner_number_of_coefficients(transformation->order); range = transformation->range; if (direction==PJ_FWD) { /* 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 */ N = *--tcy; E = *--tcx; for (; r > 0; r--) { u = *--tcy; v = *--tcx; for (c = g; 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; } static PJ_COORD horner_forward_4d (PJ_COORD point, PJ *P) { point.uv = horner_func ((HORNER *) P->opaque, PJ_FWD, point.uv); return point; } static PJ_COORD horner_reverse_4d (PJ_COORD point, PJ *P) { point.uv = horner_func ((HORNER *) P->opaque, PJ_INV, point.uv); return point; } /**********************************************************************/ static PJ_UV complex_horner (const HORNER *transformation, PJ_DIRECTION direction, PJ_UV position) { /*********************************************************************** A reimplementation of a classic Engsager/Poder Horner complex polynomial evaluation engine. ***********************************************************************/ /* 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; PJ_UV uv_error; uv_error.u = uv_error.v = HUGE_VAL; if (nullptr==transformation) return uv_error; /* 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; } /* Prepare for double Horner */ sz = 2*transformation->order + 2; range = transformation->range; 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; if (transformation->uneg) e = -e; if (transformation->vneg) n = -n; } else { /* inverse */ cb = transformation->inv_c; c = cb + sz; e = position.u - transformation->inv_origin->u; n = position.v - transformation->inv_origin->v; if (transformation->uneg) e = -e; if (transformation->vneg) n = -n; } if ((fabs(n) > range) || (fabs(e) > range)) { errno = EDOM; return uv_error; } /* 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; } position.u = E; position.v = N; return position; } static PJ_COORD complex_horner_forward_4d (PJ_COORD point, PJ *P) { point.uv = complex_horner ((HORNER *) P->opaque, PJ_FWD, point.uv); return point; } static PJ_COORD complex_horner_reverse_4d (PJ_COORD point, PJ *P) { point.uv = complex_horner ((HORNER *) P->opaque, PJ_INV, point.uv); return point; } static PJ *horner_freeup (PJ *P, int errlev) { /* Destructor */ if (nullptr==P) return nullptr; if (nullptr==P->opaque) return pj_default_destructor (P, errlev); horner_free ((HORNER *) P->opaque); P->opaque = nullptr; return pj_default_destructor (P, errlev); } static int parse_coefs (PJ *P, double *coefs, const char *param, int ncoefs) { char *buf, *init, *next = nullptr; int i; buf = static_cast<char*>(pj_calloc (strlen (param) + 2, sizeof(char))); if (nullptr==buf) { proj_log_error (P, "Horner: No memory left"); return 0; } sprintf (buf, "t%s", param); if (0==pj_param (P->ctx, P->params, buf).i) { pj_dealloc (buf); 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) { if ( next == nullptr || ','!=*next) { proj_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, complex_polynomia = 0; HORNER *Q; P->fwd4d = horner_forward_4d; P->inv4d = horner_reverse_4d; P->fwd3d = nullptr; P->inv3d = nullptr; P->fwd = nullptr; P->inv = nullptr; P->left = P->right = PJ_IO_UNITS_PROJECTED; P->destructor = horner_freeup; /* Polynomial degree specified? */ if (pj_param (P->ctx, P->params, "tdeg").i) { /* degree specified? */ degree = pj_param(P->ctx, P->params, "ideg").i; if (degree < 0 || degree > 10000) { /* What are reasonable minimum and maximums for degree? */ proj_log_debug (P, "Horner: Degree is unreasonable: %d", degree); return horner_freeup (P, PJD_ERR_INVALID_ARG); } } else { proj_log_debug (P, "Horner: Must specify polynomial degree, (+deg=n)"); return horner_freeup (P, PJD_ERR_MISSING_ARGS); } if (pj_param (P->ctx, P->params, "tfwd_c").i || pj_param (P->ctx, P->params, "tinv_c").i) /* complex polynomium? */ complex_polynomia = 1; Q = horner_alloc (degree, complex_polynomia); if (Q == nullptr) return horner_freeup (P, ENOMEM); P->opaque = Q; if (complex_polynomia) { /* Westings and/or southings? */ Q->uneg = pj_param_exists (P->params, "uneg") ? 1 : 0; Q->vneg = pj_param_exists (P->params, "vneg") ? 1 : 0; n = 2*degree + 2; if (0==parse_coefs (P, Q->fwd_c, "fwd_c", n)) return horner_freeup (P, PJD_ERR_MISSING_ARGS); if (0==parse_coefs (P, Q->inv_c, "inv_c", n)) return horner_freeup (P, PJD_ERR_MISSING_ARGS); P->fwd4d = complex_horner_forward_4d; P->inv4d = complex_horner_reverse_4d; } else { n = horner_number_of_coefficients (degree); if (0==parse_coefs (P, Q->fwd_u, "fwd_u", n)) return horner_freeup (P, PJD_ERR_MISSING_ARGS); if (0==parse_coefs (P, Q->fwd_v, "fwd_v", n)) return horner_freeup (P, PJD_ERR_MISSING_ARGS); if (0==parse_coefs (P, Q->inv_u, "inv_u", n)) return horner_freeup (P, PJD_ERR_MISSING_ARGS); if (0==parse_coefs (P, Q->inv_v, "inv_v", n)) return horner_freeup (P, PJD_ERR_MISSING_ARGS); } if (0==parse_coefs (P, (double *)(Q->fwd_origin), "fwd_origin", 2)) return horner_freeup (P, PJD_ERR_MISSING_ARGS); if (0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2)) return horner_freeup (P, PJD_ERR_MISSING_ARGS); if (0==parse_coefs (P, &Q->range, "range", 1)) Q->range = 500000; return P; }