HeteQuadrature/algoim/tanhsinh.hpp

#ifndef ALGOIM_TANHSINH_HPP
#define ALGOIM_TANHSINH_HPP

// algoim::TanhSinhQuadrature implements basic tanh-sinh quadrature methods. Some
// of the design choices follow the discussion in the paper
//    R. I. Saye, High-order quadrature on multi-component domains implicitly defined
//    by multivariate polynomials, Journal of Computational Physics, 448, 110720 (2022),
//    https://doi.org/10.1016/j.jcp.2021.110720

#include <cassert>
#include <array>
#include <vector>
#include "real.hpp"
#include "uvector.hpp"

namespace algoim
{
    struct TanhSinhQuadrature
    {
        static constexpr int n_max = 100;

        // Computes and caches all quadrature points and nodes for schemes
        // containing up to n_max many quadrature points
        static const std::array<real,n_max*(n_max+1)>& data()
        {
            auto generator = []()
            {
                std::array<real,n_max*(n_max+1)> d;
                for (int n = 1; ; ++n)
                {
                    // Generate a scheme using at most n points
                    std::vector<real> scheme(2*n);
                    int m = generate(n, &scheme[0]);
                    // m is the number of effective nodes; store the scheme
                    // which has the largest n for an effective m
                    if (m > n_max + 10)
                        break;
                    if (m <= n_max)
                        for (int i = 0; i < 2*m; ++i)
                            d[m*(m-1) + i] = scheme[i];
                }
                return d;
            };
            static auto data_ = generator();
            return data_;
        };

        // Quadrature node, relative to the interval [-1,1], for an n-point scheme
        static real x(int n, int i)
        {
            assert(1 <= n && n <= n_max && 0 <= i && i < n);
            return data()[n*(n-1) + 2*i + 0];
        }

        // Quadrature weight, relative to the interval [-1,1], for an n-point scheme
        static real w(int n, int i)
        {
            assert(1 <= n && n <= n_max && 0 <= i && i < n);
            return data()[n*(n-1) + 2*i + 1];
        }

        // Quadrature node, relative to the interval [a,b], for an n-point scheme
        static real x(int n, int i, real a, real b)
        {
            return (a + b + (b - a) * x(n, i)) * 0.5;
        }

        // Quadrature weight, relative to the interval [a,b], for an n-point scheme
        static real w(int n, int i, real a, real b)
        {
            return w(n, i) * (b - a) * 0.5;
        }

        // Generate quadrature points and weights for an n-point scheme; the result is stored
        // in pairs (x,w) in the output buffer, required to be of length at least 2n. Depending
        // on a tolerance, one or more endpoint nodes may be snapped to -1 or 1, and if more
        // than one, their quadrature weights accumulated; the result returned is the
        // reduced number of effective nodes
        static int generate(int n, real* out)
        {
            assert(n >= 1);

            // midpoint rule for one quadrature point
            if (n == 1)
            {
                out[0] = 0.0;
                out[1] = 2;
                return 1;
            }

            // Newton's method to compute the Lambert W function W(z), accurate to full precision
            auto LambertW = [](real z)
            {
                using std::log;
                using std::exp;
                real w = z < 1.0 ? z - 0.45 * z * z : 0.75 * log(z);
                // 6 iterations needed for double; 8 generally enough for quad-double;
                // definitely possible to improve with a more sophisticated initial guess
                for (int i = 0; i < 10; ++i)
                    w -= (w * exp(w) - z) / (exp(w) + w * exp(w));
                return w;
            };

            real h_a = 0.6 * 0.5 * util::pi;
            real h = 2.0 * LambertW(2.0 * h_a * (n - 1)) / n;

            int count = 0;

            // Generate central node when number of points is odd
            if (n % 2)
            {
                out[count++] = 0.0;
                out[count++] = util::pi * 0.5;
            }

            bool snappedEndpoint = false;

            // Generate remaining points symmetrically in pairs
            using std::exp;
            using std::abs;
            for (int i = 0; i < n/2; ++i)
            {
                real t = (n % 2) ? (i + 1) * h : (i + 0.5) * h;
                real exp_t = exp(t);
                real exp_mt = 1.0 / exp_t;
                real omega = 0.25 * util::pi * (exp_t - exp_mt); // omega = 0.5 pi sinh(t)
                real exp_omega = exp(omega);
                real cosh_omega_2 = exp_omega + 1.0 / exp_omega; // 2 cosh(omega)
                real cosh_t_2 = exp_t + exp_mt; // 2 cosh(t)

                // The quadrature weight function w(t) satisfies
                //   w(t) = pi (2 cosh(t)) / (2 cosh(omega))^2
                // The complementary quadrature node, 1 - x, satisfies
                //   1 - x(t) = 2 / (1 + exp(2 omega))
                real w = util::pi * cosh_t_2 / (cosh_omega_2 * cosh_omega_2);
                real xc = 2.0 / (1.0 + exp_omega * exp_omega);

                // If quadrature node is sufficiently close to -1 or 1, declare it to be exactly
                // -1 or 1 and accumulate the quad weights
                volatile real test = real(1) - xc;
                if (abs(test - real(1)) <= real(0.0))
                {
                    if (snappedEndpoint)
                    {
                        out[count-3] += w;
                        out[count-1] += w;
                    }
                    else
                    {
                        out[count++] = -1.0;
                        out[count++] = w;
                        out[count++] = 1.0;
                        out[count++] = w;
                        snappedEndpoint = true;
                    }
                }
                else
                {
                    assert(!snappedEndpoint);
                    out[count++] = -1.0 + xc;
                    out[count++] = w;
                    out[count++] = 1.0 - xc;
                    out[count++] = w;
                }
            }

            assert(count % 2 == 0 && (snappedEndpoint && count <= 2*n || !snappedEndpoint && count == 2*n));

            // Normalise quadrature weights
            real sum = 0.0;
            for (int i = 0; i < count/2; ++i)
                sum += out[2*i + 1];
            real alpha = 2.0 / sum;
            for (int i = 0; i < count/2; ++i)
                out[2*i + 1] *= alpha;

            return count/2;
        }
    };
} // namespace algoim

#endif