Chebyshev interpolation on erfc

chebyshev.cpp

#include <array>
#include <cmath>
#include <iostream>
#include <ccomplex>
#include <fftw3.h>

namespace {

template<typename>
class Plan;

template<>
class Plan<long double>
{
  private:
    ::fftwl_plan _handle;

  public:
    Plan(::fftwl_plan handle) noexcept
      : _handle(handle)
    {}

    ~Plan()
    {
        ::fftwl_destroy_plan(_handle);
    }

    void operator()() const
    {
        ::fftwl_execute(_handle);
    }
};

} // namespace

template<typename Scalar>
static Scalar chebyshev(int n, Scalar x)
{
    return std::cos(n * std::acos(x));
}

template<typename Scalar>
static Scalar target(Scalar t)
{
    Scalar x = 2 / t - 2;

    return x * x + std::log(std::erfc(x) / t);
}

template<typename Scalar>
static Scalar shift(Scalar t)
{
    return Scalar(1) / 7 * (16 * t - 9);
}

template<typename Scalar>
static Scalar unshift(Scalar cos)
{
    return (7 * cos + 9) / 16;
}

template<typename Scalar, int degree>
static std::array<Scalar, degree + 1> run()
{
    const Scalar pi = M_PIl;
    std::array<Scalar, degree + 1> c;

    for (int k = 0; k <= degree; ++k) {
        Scalar angle = pi / (degree + 1) * (k + Scalar(0.5));
        c[k] = Scalar(1) / (degree + 1) * target(unshift(std::cos(angle)));
    }

    Plan<long double>(::fftwl_plan_r2r_1d(degree + 1, &c.front(), &c.front(), FFTW_REDFT10, FFTW_ESTIMATE))();

    c[0] /= 2;

    for (int k = 0; k <= degree; ++k)
        std::cout << std::showpos << c[k] << " * chebyshev_t("
            << std::noshowpos << k << ", 2*x + 1)\n";

    return c;
}

template<int degree, typename Scalar>
static Scalar poly(Scalar t)
{
    static auto c = run<Scalar, degree>();
    Scalar r = -0.0;

    for (int k = 0; k <= degree; ++k)
        r += c[k] * chebyshev(k, shift(t));

    return r;
}

template<int degree, typename Scalar>
static Scalar eval(Scalar x)
{
    Scalar t = 1 / (1 + 0.5 * std::abs(x));
    return t * std::exp(poly<degree>(t) - x * x);
}

#include <cfloat>

int main()
{
    std::ios_base::sync_with_stdio(false);
    std::cout.precision(20);

    const int repeats = 123456;
    long double max = -0.0;

    for (int k = 0; k < repeats; ++k) {
        long double x = 16 * k / static_cast<long double>(repeats);
        long double diff = std::abs((eval<9>(x) - std::erfc(x)) / std::erfc(x));

        max = std::max(max, diff);
    }

    std::cout << "\nPoly(+1): " << poly<9>(+1.0L)
              << "\nPoly(-1): " << poly<9>(-1.0L)
              << "\nMax error: " << std::noshowpos << max << std::endl;
}