From 32e44f95d00f259cbb00abfa3eda7e271cc14842 Mon Sep 17 00:00:00 2001 From: Jake Bowhay Date: Thu, 6 Aug 2026 16:17:50 +0100 Subject: [PATCH 1/2] MAINT: improved polynomial evaluation interface --- include/xsf/evalpoly.h | 7 ++++++- include/xsf/lambertw.h | 4 ++-- include/xsf/loggamma.h | 4 ++-- tests/xsf_tests/test_orthogonal_eval.cpp | 15 +++------------ 4 files changed, 13 insertions(+), 17 deletions(-) diff --git a/include/xsf/evalpoly.h b/include/xsf/evalpoly.h index b126fb608f..bc3a6a33fc 100644 --- a/include/xsf/evalpoly.h +++ b/include/xsf/evalpoly.h @@ -19,11 +19,12 @@ #pragma once +#include "cephes/polevl.h" #include "config.h" namespace xsf { -XSF_HOST_DEVICE inline std::complex cevalpoly(const double *coeffs, int degree, std::complex z) { +XSF_HOST_DEVICE inline std::complex evalpoly(const double *coeffs, int degree, std::complex z) { /* Evaluate a polynomial with real coefficients at a complex point. * * Uses equation (3) in section 4.6.4 of [1]. Note that it is more @@ -44,4 +45,8 @@ XSF_HOST_DEVICE inline std::complex cevalpoly(const double *coeffs, int return z * a + b; } +XSF_HOST_DEVICE inline double evalpoly(const double *coeffs, int degree, double x) { + return cephes::polevl(x, coeffs, degree); +} + } // namespace xsf diff --git a/include/xsf/lambertw.h b/include/xsf/lambertw.h index 964b60caae..ec2c1fa16a 100644 --- a/include/xsf/lambertw.h +++ b/include/xsf/lambertw.h @@ -51,7 +51,7 @@ namespace detail { double coeffs[] = {-1.0 / 3.0, 1.0, -1.0}; std::complex p = std::sqrt(2.0 * (M_E * z + 1.0)); - return cevalpoly(coeffs, 2, p); + return evalpoly(coeffs, 2, p); } XSF_HOST_DEVICE inline std::complex lambertw_pade0(std::complex z) { @@ -62,7 +62,7 @@ namespace detail { /* This only gets evaluated close to 0, so we don't need a more * careful algorithm that avoids overflow in the numerator for * large z. */ - return z * cevalpoly(num, 2, z) / cevalpoly(denom, 2, z); + return z * evalpoly(num, 2, z) / evalpoly(denom, 2, z); } XSF_HOST_DEVICE inline std::complex lambertw_asy(std::complex z, long k) { diff --git a/include/xsf/loggamma.h b/include/xsf/loggamma.h index 845e99c912..9c226d9b77 100644 --- a/include/xsf/loggamma.h +++ b/include/xsf/loggamma.h @@ -50,7 +50,7 @@ namespace detail { std::complex rz = 1.0 / z; std::complex rzz = rz / z; - return (z - 0.5) * std::log(z) - z + loggamma_HLOG2PI + rz * cevalpoly(coeffs, 7, rzz); + return (z - 0.5) * std::log(z) - z + loggamma_HLOG2PI + rz * evalpoly(coeffs, 7, rzz); } XSF_HOST_DEVICE std::complex loggamma_recurrence(std::complex z) { @@ -95,7 +95,7 @@ namespace detail { 8.2246703342411321824E-1, -5.7721566490153286061E-1}; z -= 1.0; - return z * cevalpoly(coeffs, 22, z); + return z * evalpoly(coeffs, 22, z); } } // namespace detail diff --git a/tests/xsf_tests/test_orthogonal_eval.cpp b/tests/xsf_tests/test_orthogonal_eval.cpp index 133b18701b..f2631deda4 100644 --- a/tests/xsf_tests/test_orthogonal_eval.cpp +++ b/tests/xsf_tests/test_orthogonal_eval.cpp @@ -1,7 +1,6 @@ #include "../../include/xsf/config.h" #include "../testing_utils.h" -#include #include #include @@ -36,22 +35,14 @@ std::vector multiply(const std::vector &a, const std::vector &coeffs, double x) { +template +T polyval(const std::vector &coeffs, T x) { if (coeffs.size() == 1) { return coeffs[0]; } const std::vector reversed(coeffs.rbegin(), coeffs.rend()); - return xsf::cephes::polevl(x, reversed.data(), reversed.size() - 1); -} - -std::complex polyval(const std::vector &coeffs, std::complex x) { - if (coeffs.size() == 1) { - return coeffs[0]; - } - - const std::vector reversed(coeffs.rbegin(), coeffs.rend()); - return xsf::cevalpoly(reversed.data(), reversed.size() - 1, x); + return xsf::evalpoly(reversed.data(), reversed.size() - 1, x); } double sample(double a, double b, int i) { From ad228952c729e88d034c8b78f5be5006466beb46 Mon Sep 17 00:00:00 2001 From: Jake Bowhay Date: Thu, 13 Aug 2026 17:42:16 +0100 Subject: [PATCH 2/2] add degree 0 case --- include/xsf/evalpoly.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/include/xsf/evalpoly.h b/include/xsf/evalpoly.h index bc3a6a33fc..0d04c9ecb4 100644 --- a/include/xsf/evalpoly.h +++ b/include/xsf/evalpoly.h @@ -30,6 +30,9 @@ XSF_HOST_DEVICE inline std::complex evalpoly(const double *coeffs, int d * Uses equation (3) in section 4.6.4 of [1]. Note that it is more * efficient than Horner's method. */ + if (degree == 0) { + return coeffs[0]; + } double a = coeffs[0]; double b = coeffs[1]; double r = 2 * z.real(); @@ -46,6 +49,9 @@ XSF_HOST_DEVICE inline std::complex evalpoly(const double *coeffs, int d } XSF_HOST_DEVICE inline double evalpoly(const double *coeffs, int degree, double x) { + if (degree == 0) { + return coeffs[0]; + } return cephes::polevl(x, coeffs, degree); }