This file contains implementations for some special functions which are needed by the MMM family of algorithms. More...

#include <cassert>
#include <numeric>
#include <span>
#include <utility>

Include dependency graph for specfunc.hpp:

This graph shows which files directly or indirectly include this file:

Functions
double	hzeta (double order, double x)
	Hurwitz zeta function.

double	K0 (double x)
	Modified Bessel function of second kind, order 0.

double	K1 (double x)
	Modified Bessel function of second kind, order 1.

double	LPK0 (double x)
	Modified Bessel function of second kind, order 0, low precision.

double	LPK1 (double x)
	Modified Bessel function of second kind, order 1, low precision.

std::pair< double, double >	LPK01 (double x)
	Modified Bessel functions of second kind, order 0 and 1, low precision.

double	evaluateAsTaylorSeriesAt (std::span< const double > series, double x)
	Evaluate the polynomial interpreted as a Taylor series via the Horner scheme.

double	evaluateAsChebychevSeriesAt (std::span< const double > series, double x)
	Evaluate the polynomial interpreted as a Chebychev series.

Detailed Description

This file contains implementations for some special functions which are needed by the MMM family of algorithms.

These are the modified Hurwitz zeta function and the modified Bessel functions of second kind. The implementations are based on the GSL code (see specfunc.cpp for the original GSL header).

The Hurwitz zeta function is evaluated using the Euler-Maclaurin summation formula, the Bessel functions are evaluated using several different Chebychev expansions. Both achieve a precision of nearly machine precision, which is no problem for the Hurwitz zeta function, which is only used when determining the coefficients for the modified polygamma functions.

Definition in file specfunc.hpp.

Function Documentation

◆ evaluateAsChebychevSeriesAt()

double evaluateAsChebychevSeriesAt	(	std::span< const double >	series,
		double	x
	)

inline

Evaluate the polynomial interpreted as a Chebychev series.

Requires a series with at least three coefficients, i.e. no linear approximations!

Definition at line 103 of file specfunc.hpp.

Referenced by K0(), and K1().

◆ evaluateAsTaylorSeriesAt()

double evaluateAsTaylorSeriesAt	(	std::span< const double >	series,
		double	x
	)

inline

Evaluate the polynomial interpreted as a Taylor series via the Horner scheme.

Definition at line 91 of file specfunc.hpp.

Referenced by mod_psi_even(), and mod_psi_odd().

◆ hzeta()

double hzeta	(	double	order,
		double	x
	)

Hurwitz zeta function.

This function was taken from the GSL code.

Euler-Maclaurin summation formula from [28] p. 400, with several typo corrections.

Definition at line 215 of file specfunc.cpp.

References hzeta_c, and Utils::sqr().

Referenced by preparePolygammaEven(), and preparePolygammaOdd().

◆ K0()

double K0 ( double x )

Modified Bessel function of second kind, order 0.

This function was taken from the GSL code. Precise roughly up to machine precision. It is 16 times faster than std::cyl_bessel_k. If MMM1D_MACHINE_PREC is not defined, LPK0 is used instead.

Definition at line 250 of file specfunc.cpp.

References ak02_cs, ak0_cs, bi0_cs, bk0_cs, and evaluateAsChebychevSeriesAt().

Referenced by CoulombMMM1D::pair_energy(), and CoulombMMM1D::pair_force().

◆ K1()

double K1 ( double x )

Modified Bessel function of second kind, order 1.

This function was taken from the GSL code. Precise roughly up to machine precision. If MMM1D_MACHINE_PREC is not defined, LPK1 is used instead.

Definition at line 262 of file specfunc.cpp.

References ak12_cs, ak1_cs, bi1_cs, bk1_cs, and evaluateAsChebychevSeriesAt().

Referenced by far_error(), and CoulombMMM1D::pair_force().

◆ LPK0()

double LPK0 ( double x )

Modified Bessel function of second kind, order 0, low precision.

The implementation has an absolute precision of around 10^(-14), which is comparable to the relative precision sqrt implementation of current hardware in the ranges \( ]0, 8[ \) and \( ]8, 23[ \). Above 23, the precision starts to degrade, and above 27 the result drifts and slowly converges to 96% of the real value. It is 25 times faster than std::cyl_bessel_k and 1.5 times faster than K0.

Definition at line 288 of file specfunc.cpp.

References ak01_orders, ak02_cs, ak0_cs, bi0_cs, and bk0_cs.

◆ LPK01()

std::pair< double, double > LPK01 ( double x )

Modified Bessel functions of second kind, order 0 and 1, low precision.

The implementation has an absolute precision of around 10^(-14), which is comparable to the relative precision sqrt implementation of current hardware.

Definition at line 408 of file specfunc.cpp.

References ak01_orders, ak02_cs, ak0_cs, ak12_cs, ak1_cs, bi0_cs, bi1_cs, bk0_cs, and bk1_cs.

Referenced by CoulombMMM1D::pair_force().

◆ LPK1()

double LPK1 ( double x )

Modified Bessel function of second kind, order 1, low precision.

The implementation has an absolute precision of around 10^(-14), which is comparable to the relative precision sqrt implementation of current hardware in the ranges \( ]0, 8[ \) and \( ]8, 23[ \). Above 23, the precision starts to degrade, and above 27 the result drifts and slowly converges to 111% of the real value. It is 25 times faster than std::cyl_bessel_k and 1.5 times faster than K1.

Definition at line 348 of file specfunc.cpp.

References ak01_orders, ak12_cs, ak1_cs, bi1_cs, and bk1_cs.

Functions