hll/include/CubicInterpolation-internal.hpp

/*
 * Licensed to the Apache Software Foundation (ASF) under one
 * or more contributor license agreements.  See the NOTICE file
 * distributed with this work for additional information
 * regarding copyright ownership.  The ASF licenses this file
 * to you under the Apache License, Version 2.0 (the
 * "License"); you may not use this file except in compliance
 * with the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing,
 * software distributed under the License is distributed on an
 * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
 * KIND, either express or implied.  See the License for the
 * specific language governing permissions and limitations
 * under the License.
 */

#ifndef _CUBICINTERPOLATION_INTERNAL_HPP_
#define _CUBICINTERPOLATION_INTERNAL_HPP_

#include "CubicInterpolation.hpp"

#include <string>
#include <stdexcept>

namespace datasketches {

template<typename A>
static double interpolateUsingXAndYTables(const double xArr[], const double yArr[], const int offset, const double x);

template<typename A>
static double cubicInterpolate(const double x0, const double y0, const double x1, const double y1,
                               const double x2, const double y2, const double x3, const double y3, const double x);

template<typename A>
static int findStraddle(const double xArr[], const int len, const double x);

template<typename A>
static int recursiveFindStraddle(const double xArr[], const int l, const int r, const double x);

template<typename A>
static double interpolateUsingXArrAndYStride(const double xArr[], const double yStride,
                                             const int offset, const double x);

const int numEntries = 40;

//Computed for Coupon lgK = 26 ONLY. Designed for the cubic interpolator function.
const double xArrComputed[numEntries] = {
    0.0, 1.0, 20.0, 400.0,
    8000.0, 160000.0, 300000.0, 600000.0,
    900000.0, 1200000.0, 1500000.0, 1800000.0,
    2100000.0, 2400000.0, 2700000.0, 3000000.0,
    3300000.0, 3600000.0, 3900000.0, 4200000.0,
    4500000.0, 4800000.0, 5100000.0, 5400000.0,
    5700000.0, 6000000.0, 6300000.0, 6600000.0,
    6900000.0, 7200000.0, 7500000.0, 7800000.0,
    8100000.0, 8400000.0, 8700000.0, 9000000.0,
    9300000.0, 9600000.0, 9900000.0, 10200000.0
};

//Computed for Coupon lgK = 26 ONLY. Designed for the cubic interpolator function.
const double yArrComputed[numEntries] =  {
    0.0000000000000000, 1.0000000000000000, 20.0000009437402611, 400.0003963713384110,
    8000.1589294602090376, 160063.6067763759638183, 300223.7071597663452849, 600895.5933856170158833,
    902016.8065120954997838, 1203588.4983199508860707, 1505611.8245524743106216, 1808087.9449319066479802,
    2111018.0231759352609515, 2414403.2270142501220107, 2718244.7282051891088486, 3022543.7025524540804327,
    3327301.3299219091422856, 3632518.7942584538832307, 3938197.2836029687896371, 4244337.9901093561202288,
    4550942.1100616492331028, 4858010.8438911894336343, 5165545.3961938973516226, 5473546.9757476449012756,
    5782016.7955296505242586, 6090956.0727340159937739, 6400366.0287892958149314, 6710247.8893762007355690,
    7020602.8844453142955899, 7331432.2482349723577499, 7642737.2192891482263803, 7954519.0404754765331745,
    8266778.9590033423155546, 8579518.2264420464634895, 8892738.0987390466034412, 9206439.8362383283674717,
    9520624.7036988288164139, 9835293.9703129194676876, 10150448.9097250290215015, 10466090.8000503256917000
};

template<typename A>
double CubicInterpolation<A>::usingXAndYTables(const double x) {
  return usingXAndYTables(xArrComputed, yArrComputed, numEntries, x);
}

template<typename A>
double CubicInterpolation<A>::usingXAndYTables(const double xArr[], const double yArr[],
                                            const int len, const double x) {
  int offset;
  if (x < xArr[0] || x > xArr[len-1]) {
    throw std::invalid_argument("x value out of range: " + std::to_string(x));
  }

  if (x ==  xArr[len-1]) { // corner case
    return (yArr[len-1]);
  }

  offset = findStraddle<A>(xArr, len, x);
  if (offset < 0 && offset > len-2) {
    throw std::logic_error("offset must be >= 0 and <= " + std::to_string(len) + "-2");
  }

  if (offset == 0) { // corner case
    return (interpolateUsingXAndYTables<A>(xArr, yArr, (offset-0), x));
  }
  else if (offset == numEntries-2) { // corner case
    return (interpolateUsingXAndYTables<A>(xArr, yArr, (offset-2), x));
  }
  else { // main case
    return (interpolateUsingXAndYTables<A>(xArr, yArr, (offset-1), x));
  }
  throw std::logic_error("Exception should be unreachable");
}

// In C: again-two-registers cubic_interpolate_aux L1368
template<typename A>
static double interpolateUsingXAndYTables(const double xArr[], const double yArr[],
                                          const int offset, const double x) {
    return (cubicInterpolate<A>(xArr[offset+0], yArr[offset+0],
                        xArr[offset+1], yArr[offset+1],
                        xArr[offset+2], yArr[offset+2],
                        xArr[offset+3], yArr[offset+3],
                        x) );
}

template<typename A>
static inline double cubicInterpolate(const double x0, const double y0,
                                      const double x1, const double y1,
                                      const double x2, const double y2,
                                      const double x3, const double y3,
                                      const double x)
{
  double l0_numer = (x - x1) * (x - x2) * (x - x3);
  double l1_numer = (x - x0) * (x - x2) * (x - x3);
  double l2_numer = (x - x0) * (x - x1) * (x - x3);
  double l3_numer = (x - x0) * (x - x1) * (x - x2);

  double l0_denom = (x0 - x1) * (x0 - x2) * (x0 - x3);
  double l1_denom = (x1 - x0) * (x1 - x2) * (x1 - x3);
  double l2_denom = (x2 - x0) * (x2 - x1) * (x2 - x3);
  double l3_denom = (x3 - x0) * (x3 - x1) * (x3 - x2);

  double term0 = y0 * l0_numer / l0_denom;
  double term1 = y1 * l1_numer / l1_denom;
  double term2 = y2 * l2_numer / l2_denom;
  double term3 = y3 * l3_numer / l3_denom;

  return (term0 + term1 + term2 + term3);
}

/* returns j such that xArr[j] <= x and x < xArr[j+1] */
template<typename A>
static int findStraddle(const double xArr[], const int len, const double x)
{
  if ((len < 2) || (x < xArr[0]) || (x > xArr[len-1])) {
    throw std::logic_error("invariant violated during interpolation");
  }
  return(recursiveFindStraddle<A>(xArr, 0, len-1, x));
}


/* the invariant here is that xArr[l] <= x && x < xArr[r] */
template<typename A>
static int recursiveFindStraddle(const double xArr[], const int l, const int r, const double x)
{
  int m;
  if (l >= r) {
    throw std::logic_error("lower bound not less than upper bound in search");
  }
  if ((xArr[l] > x) || (x >= xArr[r])) { // the invariant
    throw std::logic_error("target value invariant violated in search");
  }

  if (l+1 == r) return (l);
  m = l + ((r-l)/2);
  if (xArr[m] <= x) return (recursiveFindStraddle<A>(xArr, m, r, x));
  else              return (recursiveFindStraddle<A>(xArr, l, m, x));
}


//Interpolate using X table and Y stride

/**
 * Cubic interpolation using interpolation X table and Y stride.
 *
 * @param xArr The x array
 * @param yStride The y stride
 * @param xArrLen the length of xArr
 * @param x The value x
 * @return cubic interpolation
 */
//In C: again-two-registers cubic_interpolate_with_x_arr_and_y_stride L1411
// Used by HllEstimators
template<typename A>
double CubicInterpolation<A>::usingXArrAndYStride(const double xArr[], const int xArrLen,
                                                  const double yStride, const double x) {
  const int xArrLenM1 = xArrLen - 1;

  if ((xArrLen < 4) || (x < xArr[0]) || (x > xArr[xArrLenM1])) {
    throw std::logic_error("impossible values during interpolaiton");
  }

  if (x ==  xArr[xArrLenM1]) { /* corner case */
    return (yStride * (xArrLenM1));
  }

  const int offset = findStraddle<A>(xArr, xArrLen, x); //uses recursion
  const int xArrLenM2 = xArrLen - 2;
  if ((offset < 0) || (offset > xArrLenM2)) {
    throw std::logic_error("invalid offset during interpolation");
  }

  if (offset == 0) { /* corner case */
    return (interpolateUsingXArrAndYStride<A>(xArr, yStride, (offset - 0), x));
  }
  else if (offset == xArrLenM2) { /* corner case */
    return (interpolateUsingXArrAndYStride<A>(xArr, yStride, (offset - 2), x));
  }
  /* main case */
  return (interpolateUsingXArrAndYStride<A>(xArr, yStride, (offset - 1), x));
}

//In C: again-two-registers cubic_interpolate_with_x_arr_and_y_stride_aux L1402
template<typename A>
static double interpolateUsingXArrAndYStride(const double xArr[], const double yStride,
                                             const int offset, const double x) {
  return cubicInterpolate<A>(
      xArr[offset + 0], yStride * (offset + 0),
      xArr[offset + 1], yStride * (offset + 1),
      xArr[offset + 2], yStride * (offset + 2),
      xArr[offset + 3], yStride * (offset + 3),
      x);
}

}

#endif // _CUBICINTERPOLATION_INTERNAL_HPP_