cubicsuperancillary.hpp

#pragma once 
#include <vector>

namespace teqp {

namespace CubicSuperAncillary {

struct Chebyshev{
public:
    const std::vector<double> coeff;
    const double xmin, xmax;
    // Evaluate the expansion with Clenshaw's method
    double y(double x) const{
        // Scale to (-1, 1)
        double xscaled = (2*x - (xmax + xmin)) / (xmax - xmin);
        int Norder = static_cast<int>(coeff.size()) - 1;
        double u_k = 0, u_kp1 = coeff[Norder], u_kp2 = 0;
        for (int k = Norder-1; k > 0; k--){ // k must be signed!
            // Do the recurrent calculation
            u_k = 2.0*xscaled*u_kp1 - u_kp2 + coeff[k];
            // Update the values
            u_kp2 = u_kp1; u_kp1 = u_k;
        }
        return coeff[0] + xscaled*u_kp1 - u_kp2;
    };
};

// https://proquest.safaribooksonline.com/9780321637413
// https://web.stanford.edu/class/archive/cs/cs107/cs107.1202/lab1/
static int midpoint_Knuth(int x, int y) {
    return (x & y) + ((x ^ y) >> 1);
};

struct SuperAncillary{
public:

    const std::vector<Chebyshev> exps;

    int get_index(double x) const{
        int iL = 0, iR = static_cast<int>(exps.size()) - 1, iM;
        while (iR - iL > 1) {
            iM = midpoint_Knuth(iL, iR);
            if (x >= exps[iM].xmin) {
                iL = iM;
            }
            else {
                iR = iM;
            }
        }
        return (x < exps[iL].xmax) ? iL : iR;
    };

    /// Evaluate the SuperAncillary
    double y(double x) const{
        // First check whether the input is possible
        if (x < exps[0].xmin) {
            throw std::invalid_argument("Ttilde (" + std::to_string(x) + ") is below the minimum of " + std::to_string(exps[0].xmin));
        }
        if (x > exps.back().xmax) {
            throw std::invalid_argument("Ttilde (" + std::to_string(x) + ") is above the maximum of " + std::to_string(exps.back().xmax));
        }
        // Bisection to find the expansion
        // we need
        auto i = get_index(x);
        // Evaluate the expansion
        return exps[i].y(x);
    }
};

const auto vdW_p = SuperAncillary{
{
  {
    {
      8.369016487003571e-13,
      1.3838182749488314e-12,
      8.006631877160537e-13,
      3.3527912014616007e-13,
      1.043509859614108e-13,
      2.4518354339464807e-14,
      4.362454910953413e-15,
      5.800290149785568e-16,
      5.5246311563342216e-17,
      3.341712352615597e-18,
      6.943507182585673e-20,
      -6.6172300100790814e-21,
      -4.799249815587726e-22,
      8.197853768324654e-24,
      1.6972506557505968e-24,
      -1.4268299756055458e-26,
      -5.664810160392086e-27,
      4.437342591868191e-31,
      -1.4685138788754898e-28
    },
    0.02962962962962963,
    0.03796296296296296
  },
  {
    {
      1.1901434626227093e-10,
      1.7361780173369773e-10,
      7.639416392753658e-11,
      2.219331080562664e-11,
      4.468985012784936e-12,
      6.345623163747038e-13,
      6.256832040916454e-14,
      3.9794258615342e-15,
      1.2126532982344795e-16,
      -2.510288446235168e-18,
      -3.1189177665897616e-19,
      -3.0550813286427195e-22,
      6.384476478410435e-22,
      2.6783405454063793e-25,
      -1.4724057241513977e-24,
      -5.563678192342752e-26,
      1.0097419586828951e-27,
      6.679443056687351e-26,
      5.4526065768876336e-27
    },
    0.03796296296296296,
    0.0462962962962963
  },
  {
    {
      3.341667610526653e-08,
      5.118389108542411e-08,
      2.4371256768154856e-08,
      7.57899054003236e-09,
      1.5635438033324503e-09,
      2.0715013642702e-10,
      1.5165995758982126e-11,
      1.6314454982752663e-13,
      -6.017223931292203e-14,
      -1.9704412673913643e-15,
      2.8051398344015177e-16,
      7.022929896522875e-18,
      -1.5199614910962848e-18,
      1.1958046052493347e-20,
      7.00911492046808e-21,
      -4.237329623030918e-22,
      -3.863191954564062e-23,
      -1.1166938269465874e-23,
      -1.2818068320304144e-23
    },
    0.0462962962962963,
    0.06296296296296297
  },
  {
    {
      8.410724206747044e-06,
      1.2481631400263709e-05,
      5.378159443205515e-06,
      1.3704527561658967e-06,
      1.9136168235190198e-07,
      8.957711663130537e-09,
      -9.462672897450059e-10,
      -6.968811539371602e-11,
      1.0243849900317855e-11,
      2.337729169280951e-13,
      -1.1494530386099005e-13,
      5.3160592525426036e-15,
      7.530730564512324e-16,
      -1.1431990843247566e-16,
      2.805445086019535e-18,
      7.988918627843267e-19,
      -1.0295338058018748e-19,
      7.210533399624742e-21,
      5.761312966431838e-22
    },
    0.06296296296296297,
    0.0962962962962963
  },
  {
    {
      0.0006802593682753608,
      0.0009003332567452364,
      0.0002880297751978097,
      4.118964226458961e-05,
      7.986981798077676e-07,
      -2.628431475282669e-07,
      1.4516582903911464e-08,
      1.5647152494426528e-09,
      -4.1929621787396065e-10,
      3.8557343395626107e-11,
      8.368753107866514e-13,
      -8.31125296022126e-13,
      1.3326079144488354e-13,
      -8.864161253017153e-15,
      -8.83270499045497e-16,
      3.3879136780332834e-16,
      -4.8188763160481214e-17,
      2.9076523496871995e-18,
      3.8719993601362204e-19
    },
    0.0962962962962963,
    0.16296296296296295
  },
  {
    {
      0.015871052250174933,
      0.017373046113870208,