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