Skip to content

MAINT: improved polynomial evaluation interface #243

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fbourgey merged 2 commits into from
Aug 13, 2026
Merged

MAINT: improved polynomial evaluation interface #243

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
32e44f9
MAINT: improved polynomial evaluation interface
j-bowhay Aug 6, 2026
ad22895
add degree 0 case
j-bowhay Aug 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 6 additions & 1 deletion include/xsf/evalpoly.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,11 +19,12 @@

#pragma once

#include "cephes/polevl.h"
#include "config.h"

namespace xsf {

XSF_HOST_DEVICE inline std::complex<double> cevalpoly(const double *coeffs, int degree, std::complex<double> z) {
XSF_HOST_DEVICE inline std::complex<double> evalpoly(const double *coeffs, int degree, std::complex<double> z) {

@fbourgey fbourgey Aug 7, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about if degree=0? Should we add a guard

if (degree=0) {
   return coeffs[0];
}

as coeffs[1] might not exist. The same would apply for the float overload.

@j-bowhay j-bowhay Aug 13, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done - there is quite a bit else that could be done such as checking degree is non-negative and adding standalone tests but I would rather not grow the scope of this PR

/* 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
Expand All @@ -44,4 +45,8 @@ XSF_HOST_DEVICE inline std::complex<double> 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
4 changes: 2 additions & 2 deletions include/xsf/lambertw.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ namespace detail {
double coeffs[] = {-1.0 / 3.0, 1.0, -1.0};
std::complex<double> 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<double> lambertw_pade0(std::complex<double> z) {
Expand All @@ -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<double> lambertw_asy(std::complex<double> z, long k) {
Expand Down
4 changes: 2 additions & 2 deletions include/xsf/loggamma.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ namespace detail {
std::complex<double> rz = 1.0 / z;
std::complex<double> 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<double> loggamma_recurrence(std::complex<double> z) {
Expand Down Expand Up @@ -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

Expand Down
15 changes: 3 additions & 12 deletions tests/xsf_tests/test_orthogonal_eval.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
#include "../../include/xsf/config.h"
#include "../testing_utils.h"

#include <xsf/cephes/polevl.h>
#include <xsf/evalpoly.h>
#include <xsf/orthogonal_eval.h>

Expand Down Expand Up @@ -36,22 +35,14 @@ std::vector<double> multiply(const std::vector<double> &a, const std::vector<dou
return out;
}

double polyval(const std::vector<double> &coeffs, double x) {
template <typename T>
T polyval(const std::vector<double> &coeffs, T x) {
if (coeffs.size() == 1) {
return coeffs[0];
}

const std::vector<double> reversed(coeffs.rbegin(), coeffs.rend());
return xsf::cephes::polevl(x, reversed.data(), reversed.size() - 1);
}

std::complex<double> polyval(const std::vector<double> &coeffs, std::complex<double> x) {
if (coeffs.size() == 1) {
return coeffs[0];
}

const std::vector<double> 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) {
Expand Down
Loading