Pol

CMakeLists.txt

cmake_minimum_required(VERSION 3.25)
project(polyinterp)

add_executable(polyinterp polyinterp.cpp)

int.c

/* polint.f -- translated by f2c (version 19970805).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* DECK POLINT */
/* Subroutine */ int polint_(integer *n, real *x, real *y, real *c__, integer 
	*ierr)
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    integer i__, k, km1;
    real dif;

/* ***BEGIN PROLOGUE  POLINT */
/* ***PURPOSE  Produce the polynomial which interpolates a set of discrete
 */
/*            data points. */
/* ***LIBRARY   SLATEC */
/* ***CATEGORY  E1B */
/* ***TYPE      SINGLE PRECISION (POLINT-S, DPLINT-D) */
/* ***KEYWORDS  POLYNOMIAL INTERPOLATION */
/* ***AUTHOR  Huddleston, R. E., (SNLL) */
/* ***DESCRIPTION */

/*     Written by Robert E. Huddleston, Sandia Laboratories, Livermore */

/*     Abstract */
/*        Subroutine POLINT is designed to produce the polynomial which */
/*     interpolates the data  (X(I),Y(I)), I=1,...,N.  POLINT sets up */
/*     information in the array C which can be used by subroutine POLYVL 
*/
/*     to evaluate the polynomial and its derivatives and by subroutine */
/*     POLCOF to produce the coefficients. */

/*     Formal Parameters */
/*     N  - the number of data points  (N .GE. 1) */
/*     X  - the array of abscissas (all of which must be distinct) */
/*     Y  - the array of ordinates */
/*     C  - an array of information used by subroutines */
/*     *******  Dimensioning Information  ******* */
/*     Arrays X,Y, and C must be dimensioned at least N in the calling */
/*     program. */

/* ***REFERENCES  L. F. Shampine, S. M. Davenport and R. E. Huddleston, */
/*                 Curve fitting by polynomials in one variable, Report */
/*                 SLA-74-0270, Sandia Laboratories, June 1974. */
/* ***ROUTINES CALLED  XERMSG */
/* ***REVISION HISTORY  (YYMMDD) */
/*   740601  DATE WRITTEN */
/*   861211  REVISION DATE from Version 3.2 */
/*   891214  Prologue converted to Version 4.0 format.  (BAB) */
/*   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ) */
/*   920501  Reformatted the REFERENCES section.  (WRB) */
/* ***END PROLOGUE  POLINT */
/* ***FIRST EXECUTABLE STATEMENT  POLINT */
    /* Parameter adjustments */
    --c__;
    --y;
    --x;

    /* Function Body */
    if (*n <= 0) {
	goto L91;
    }
    c__[1] = y[1];
    if (*n == 1) {
	return 0;
    }
    i__1 = *n;
    for (k = 2; k <= i__1; ++k) {
	c__[k] = y[k];
	km1 = k - 1;
	i__2 = km1;
	for (i__ = 1; i__ <= i__2; ++i__) {
/*     CHECK FOR DISTINCT X VALUES */
	    dif = x[i__] - x[k];
	    if (dif == 0.f) {
		goto L92;
	    }
	    c__[k] = (c__[i__] - c__[k]) / dif;
/* L10010: */
	}
    }
    *ierr = 0;
    return 0;
L91:
    *ierr = 1;
/*     N IS ZERO OR NEGATIVE */
    return 0;
L92:
    *ierr = 2;
/*     THE ABSCISSAS ARE NOT DISTINCT */
    return 0;
} /* polint_ */

polyinterp.cpp

#define TEST
#include "polyinterp.h"

polyinterp.h

// -*- c++ -*- $Id: polyinterp.h,v 1.2 2000/11/01 11:19:05 royr Exp $
// SPDX-License-Identifier: MIT
//
// A polynomial interpolator.
//
// Operation add(x,y) adds point (x,y) to the interpolating polynomial.
// All the abscissae must be unique.  No two abscissae can be alike.
// How close can they be?  A threshold specifies the minimum distance.
// Points with abscissae closer than the minimum threshold merge
// at the arithmetic mean.

#ifdef __cplusplus
extern "C" {
#endif
#include "slatec_polint.h"
#include "slatec_polyvl.h"
#ifdef __cplusplus
}
#endif

#ifdef __cplusplus

#include <functional>
#include <limits>
#include <vector>

template <typename Scalar>
enum slatec_polint_status polint(size_t n, const Scalar x[], const Scalar y[],
                                 Scalar c[]);

template <typename Scalar>
enum slatec_polyvl_status polyvl(Scalar xx, Scalar *yy, size_t n,
                                 const Scalar x[], const Scalar c[]);

template <typename Scalar>
struct poly_interpolator // a unary functor
{
private:
  Scalar abscissaDeltaThres;
  std::vector<Scalar> X, Y, C;
  std::vector<int> N;

public:
  poly_interpolator() : abscissaDeltaThres(0) {}

  void set_abscissa_thres(Scalar const &x) {
    if (0 <= x)
      abscissaDeltaThres = x;
  }

  void add(Scalar const &x, Scalar const &y) {
    // sort x by insertion -- iterate while X[i] < x
    auto Xi = X.begin();
    for (; Xi != X.end() && *Xi < x; ++Xi)
      ;
    // Xi == X.end() || *Xi >= x
    auto i = std::distance(X.begin(), Xi);
    if (Xi != X.begin() && x - Xi[-1] <= abscissaDeltaThres) {
      --i;
      X[i] = (x + X[i] * N[i]) / (N[i] + 1);
      Y[i] = (y + Y[i] * N[i]) / (N[i] + 1);
      ++N[i];
    } else if (Xi != X.end() && Xi[0] - x <= abscissaDeltaThres) {
      X[i] = (x + X[i] * N[i]) / (N[i] + 1);
      Y[i] = (y + Y[i] * N[i]) / (N[i] + 1);
      ++N[i];
    } else {
      // Get the dangerous bit over with!  Throwing is the worry.
      // It's fine---except half way through an add op.  Since there
      // are four separate vectors to update.  Four memory
      // re-allocations aren't atomic!
      X.reserve(N.size() + 1);
      Y.reserve(N.size() + 1);
      C.reserve(N.size() + 1);
      N.reserve(N.size() + 1);

      Xi = X.begin();
      auto Yi = Y.begin();
      auto Ci = C.begin();
      auto Ni = N.begin();
      std::advance(Xi, i);
      std::advance(Yi, i);
      std::advance(Ci, i);
      std::advance(Ni, i);
      X.insert(Xi, x);
      Y.insert(Yi, y);
      C.insert(Ci, 0);
      N.insert(Ni, 1);
    }
  }

  void interpolate() {
    enum slatec_polint_status status =
        polint(N.size(), X.data(), Y.data(), C.data());
    if (status != slatec_polint_success)
      throw status;
  }

  Scalar operator()(const Scalar &x) const {
    if (N.size() == 0)
      return x;
    Scalar y;
    enum slatec_polyvl_status status =
        polyvl(x, &y, N.size(), X.data(), C.data());
    if (status != slatec_polyvl_success)
      throw status;
    return y;
  }

  size_t n() const   // How many interpolating
  {                  // points evaluate the
    return N.size(); // polynomial?
  }

  void clear() {
    X.clear();
    Y.clear();
    C.clear();
    N.clear();
  }
};

////////////////////////////////////////////////////////////////////////

#ifdef TEST

#include <stdlib.h>
#include <unistd.h>

#include <cstdio>

int main(int argc, char *argv[]) {
  poly_interpolator<float> poly;
  double a = -1, b = 1, c = 0.1;
  int opt;
  while ((opt = getopt(argc, argv, "a:b:c:d:")) != -1)
    switch (opt) {
    case 'a':
      a = atof(optarg);
      break;
    case 'b':
      b = atof(optarg);
      break;
    case 'c':
      c = atof(optarg);
      break;
    case 'd':
      poly.set_abscissa_thres(atof(optarg));
    }
  double x, y;
  while (optind < argc && sscanf(argv[optind++], " %lf,%lf", &x, &y) == 2)
    poly.add(x, y);
  poly.interpolate();
  for (x = a; x < b; x += c)
    printf("%lf,%lf\n", x, poly(x));
  return EXIT_SUCCESS;
}

// g++ -o polyinterp -O2 -DTEST -x c++ polyinterp.h

#endif

template <>
enum slatec_polint_status polint<double>(size_t n, const double x[],
                                         const double y[], double c[]) {
  return slatec_polint(n, x, y, c);
}

template <>
enum slatec_polyvl_status polyvl<double>(double xx, double *yy, size_t n,
                                         const double x[], const double c[]) {
  return slatec_polyvl(xx, yy, n, x, c);
}

template <>
enum slatec_polint_status polint<float>(size_t n, const float x[],
                                        const float y[], float c[]) {
  return slatec_polintf(n, x, y, c);
}

template <>
enum slatec_polyvl_status polyvl<float>(float xx, float *yy, size_t n,
                                        const float x[], const float c[]) {
  return slatec_polyvlf(xx, yy, n, x, c);
}

#endif // __cplusplus

slatec_polint.h

#pragma once

/*
 * for size_t
 */
#include <stddef.h>

enum slatec_polint_status {
  slatec_polint_success,
  slatec_polint_failure,
  slatec_polint_abscissae_not_distinct
};

/*!
 * \brief Double-precision polynomial interpolation.
 *
 * \details Function \c slatec_polint generates a polynomial that interpolates
 * the vectors \c x[i] by \c y[i] where \c i ranges from 1 to \c{n}. It prepares
 * information in the array \c c that can be utilized by the subroutine \c
 * polyvl to calculate the polynomial and its derivatives.
 *
 * \see https://netlib.org/slatec/src/polint.f
 */
static inline enum slatec_polint_status
slatec_polint(size_t n, const double x[], const double y[], double c[]) {
  if (n == 0)
    return slatec_polint_failure;
  c[0] = y[0];
  if (n == 1)
    return slatec_polint_success;
  for (size_t k = 1; k < n; k++) {
    c[k] = y[k];
    for (size_t i = 0; i < k; i++) {
      const double dif = x[i] - x[k];
      if (dif == 0)
        return slatec_polint_abscissae_not_distinct;
      c[k] = (c[i] - c[k]) / dif;
    }
  }
  return slatec_polint_success;
}

/*!
 * \brief Single-precision polynomial interpolation.
 */
static inline enum slatec_polint_status
slatec_polintf(size_t n, const float x[], const float y[], float c[]) {
  if (n == 0)
    return slatec_polint_failure;
  c[0] = y[0];
  if (n == 1)
    return slatec_polint_success;
  for (size_t k = 1; k < n; k++) {
    c[k] = y[k];
    for (size_t i = 0; i < k; i++) {
      const float dif = x[i] - x[k];
      if (dif == 0)
        return slatec_polint_abscissae_not_distinct;
      c[k] = (c[i] - c[k]) / dif;
    }
  }
  return slatec_polint_success;
}

slatec_polyvl.h

#pragma once

/*
 * for size_t
 */
#include <stddef.h>

enum slatec_polyvl_status { slatec_polyvl_success, slatec_polyvl_failure };

/*!
 * \brief Computes polynomial double-precision values.
 *
 * \details Calculates the value of a polynomial where the polynomial was
 * produced by a previous call to \c{polint}. The argument \c n and the arrays
 * \c x and \c c must not be altered between the call to \c polint and the call
 * to \c{slatec_polyvl}.
 *
 * Simplifies the original Fortran.
 *
 * \see https://netlib.org/slatec/src/polyvl.f
 */
static inline enum slatec_polyvl_status slatec_polyvl(double xx, double *yy,
                                                      size_t n,
                                                      const double x[],
                                                      const double c[]) {
  if (n == 0)
    return slatec_polyvl_failure;
  double pione = 1, pone = c[0];
  if (n == 1) {
    *yy = pone;
    return slatec_polyvl_success;
  }
  double ptwo;
  /*
   * This makes an assumption about the number of iterations for the k=[1, n)
   * loop. For \c ptwo to receive a value, at least one iteration must run.
   * The compiler may give warnings.
   */
  for (size_t k = 1; k < n; k++) {
    const double pitwo = (xx - x[k - 1]) * pione;
    pione = pitwo;
    ptwo = pone + pitwo * c[k];
    pone = ptwo;
  }
  *yy = ptwo;
  return slatec_polyvl_success;
}

/*!
 * \brief Computes polynomial single-precision values.
 */
static inline enum slatec_polyvl_status slatec_polyvlf(float xx, float *yy,
                                                       size_t n,
                                                       const float x[],
                                                       const float c[]) {
  if (n == 0)
    return slatec_polyvl_failure;
  float pione = 1, pone = c[0];
  if (n == 1) {
    *yy = pone;
    return slatec_polyvl_success;
  }
  float ptwo;
  for (size_t k = 1; k < n; k++) {
    const float pitwo = (xx - x[k - 1]) * pione;
    pione = pitwo;
    ptwo = pone + pitwo * c[k];
    pone = ptwo;
  }
  *yy = ptwo;
  return slatec_polyvl_success;
}

yvl.c

/* polyvl.f -- translated by f2c (version 19970805).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* DECK POLYVL */
/* Subroutine */ int polyvl_(real *xx, real *yfit, integer *n, real *x, real *
	c__, integer *ierr)
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    real pone, ptwo;
    integer k;
    real pione, pitwo;

/* ***BEGIN PROLOGUE  POLYVL */
/* ***PURPOSE  Calculate the value of a polynomial and its first NDER */
/*            derivatives where the polynomial was produced by a previous 
*/
/*            call to POLINT. */
/* ***LIBRARY   SLATEC */
/* ***CATEGORY  E3 */
/* ***TYPE      SINGLE PRECISION (POLYVL-S, DPOLVL-D) */
/* ***KEYWORDS  POLYNOMIAL EVALUATION */
/* ***AUTHOR  Huddleston, R. E., (SNLL) */
/* ***DESCRIPTION */

/*     Written by Robert E. Huddleston, Sandia Laboratories, Livermore */

/*     Abstract - */
/*        Subroutine POLYVL calculates the value of the polynomial and */
/*     its first NDER derivatives where the polynomial was produced by */
/*     a previous call to POLINT. */
/*        The variable N and the arrays X and C must not be altered */
/*     between the call to POLINT and the call to POLYVL. */

/*     ******  Dimensioning Information ******* */

/*     YP   must be dimensioned by at least NDER */
/*     X    must be dimensioned by at least N (see the abstract ) */
/*     C    must be dimensioned by at least N (see the abstract ) */
/*     WORK must be dimensioned by at least 2*N if NDER is .GT. 0. */

/*     *** Note *** */
/*       If NDER=0, neither YP nor WORK need to be dimensioned variables. 
*/
/*       If NDER=1, YP does not need to be a dimensioned variable. */


/*     *****  Input parameters */

/*     NDER - the number of derivatives to be evaluated */

/*     XX   - the argument at which the polynomial and its derivatives */
/*            are to be evaluated. */

/*     N    - ***** */
/*            *       N, X, and C must not be altered between the call */
/*     X    - *       to POLINT and the call to POLYVL. */
/*     C    - ***** */


/*     *****  Output Parameters */

/*     YFIT - the value of the polynomial at XX */

/*     YP   - the derivatives of the polynomial at XX.  The derivative of 
*/
/*            order J at XX is stored in  YP(J) , J = 1,...,NDER. */

/*     IERR - Output error flag with the following possible values. */
/*          = 1  indicates normal execution */

/*     ***** Storage Parameters */

/*     WORK  = this is an array to provide internal working storage for */
/*             POLYVL.  It must be dimensioned by at least 2*N if NDER is 
*/
/*             .GT. 0.  If NDER=0, WORK does not need to be a dimensioned 
*/
/*             variable. */

/* ***REFERENCES  L. F. Shampine, S. M. Davenport and R. E. Huddleston, */
/*                 Curve fitting by polynomials in one variable, Report */
/*                 SLA-74-0270, Sandia Laboratories, June 1974. */
/* ***ROUTINES CALLED  (NONE) */
/* ***REVISION HISTORY  (YYMMDD) */
/*   740601  DATE WRITTEN */
/*   890531  Changed all specific intrinsics to generic.  (WRB) */
/*   890531  REVISION DATE from Version 3.2 */
/*   891214  Prologue converted to Version 4.0 format.  (BAB) */
/*   920501  Reformatted the REFERENCES section.  (WRB) */
/* ***END PROLOGUE  POLYVL */
/* ***FIRST EXECUTABLE STATEMENT  POLYVL */
    /* Parameter adjustments */
    --c__;
    --x;

    /* Function Body */
    *ierr = 1;

/*     *****   CODING FOR THE CASE NDER = 0 */

    pione = 1.f;
    pone = c__[1];
    *yfit = pone;
    if (*n == 1) {
	return 0;
    }
    i__1 = *n;
    for (k = 2; k <= i__1; ++k) {
	pitwo = (*xx - x[k - 1]) * pione;
	pione = pitwo;
	ptwo = pone + pitwo * c__[k];
	pone = ptwo;
/* L10010: */
    }
    *yfit = ptwo;
    return 0;

/*     *****   END OF NDER = 0 CASE */

} /* polyvl_ */