Implementation of specfunc.hpp. More...

#include "specfunc.hpp"
#include <utils/math/sqr.hpp>
#include <cmath>
#include <numbers>
#include <tuple>
#include <utility>

Include dependency graph for specfunc.cpp:

Go to the source code of this file.

Functions
double	hzeta (double s, double q)
	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.

Variables
static double const	hzeta_c [15]
	Coefficients for Maclaurin summation in hzeta().

static int	ak01_orders []
	necessary orders for K0/1 from 2 up to 22 for 10^-14 precision.

Chebyshev expansions based on SLATEC bk0(), bk0e()
static double	bk0_cs [11]
	Series for `bk0`.

static double	ak0_cs [17]
	Series for `ak0`.

static double	ak02_cs [14]
	Series for `ak02`.

Chebyshev expansions based on SLATEC besi0()
static double	bi0_cs [12]
	Series for `bi0`.

Chebyshev expansions based on SLATEC besk1(), besk1e()
static double	bk1_cs [11]
	Series for `bk1`.

static double	ak1_cs [17]
	Series for `ak1`.

static double	ak12_cs [14]
	Series for `ak12`.

Chebyshev expansions based on SLATEC besi1(), besi1e()
static double	bi1_cs [11]
	Series for `bi1`.

Detailed Description

Implementation of specfunc.hpp.

Definition in file specfunc.cpp.

Function Documentation

◆ 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.

Variable Documentation

◆ ak01_orders

int ak01_orders[]

static

Initial value:

= {
    
    11, 11, 10, 10, 9, 9,
    
    6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 2, 2, 2, 2, 2}

necessary orders for K0/1 from 2 up to 22 for 10^-14 precision.

Note that at 8 the expansion changes. From 23 to 26 order 2 is used, above order 1. For the latter cases, separate implementations are necessary.

Definition at line 282 of file specfunc.cpp.

Referenced by LPK0(), LPK01(), and LPK1().

◆ ak02_cs

double ak02_cs[14]

static

Initial value:

= {
    2.5 - 0.01201869826307592, -0.00917485269102569, 0.00014445509317750,
    -0.00000401361417543,      0.00000015678318108,  -0.00000000777011043,
    0.00000000046111825,       -0.00000000003158592, 0.00000000000243501,
    -0.00000000000020743,      0.00000000000001925,  -0.00000000000000192,
    0.00000000000000020,       -0.00000000000000002}

Series for ak02.

On the interval 0. to 1.25000d-01

Label	Value
with weighted error	2.34e-17
log weighted error	16.63
significant figures required	14.67
decimal places required	17.20

Definition at line 107 of file specfunc.cpp.

Referenced by K0(), LPK0(), and LPK01().

◆ ak0_cs

double ak0_cs[17]

static

Initial value:

= {
    2.5 - 0.07643947903327941, -0.02235652605699819, 0.00077341811546938,
    -0.00004281006688886,      0.00000308170017386,  -0.00000026393672220,
    0.00000002563713036,       -0.00000000274270554, 0.00000000031694296,
    -0.00000000003902353,      0.00000000000506804,  -0.00000000000068895,
    0.00000000000009744,       -0.00000000000001427, 0.00000000000000215,
    -0.00000000000000033,      0.00000000000000005}

Series for ak0.

On the interval 1.25000d-01 to 5.00000d-01

Label	Value
with weighted error	5.34e-17
log weighted error	16.27
significant figures required	14.92
decimal places required	16.89

Definition at line 90 of file specfunc.cpp.

Referenced by K0(), LPK0(), and LPK01().

◆ ak12_cs

double ak12_cs[14]

static

Initial value:

= {
    2.5 + 0.06379308343739001, 0.02832887813049721,  -0.00024753706739052,
    0.00000577197245160,       -0.00000020689392195, 0.00000000973998344,
    -0.00000000055853361,      0.00000000003732996,  -0.00000000000282505,
    0.00000000000023720,       -0.00000000000002176, 0.00000000000000215,
    -0.00000000000000022,      0.00000000000000002}

Series for ak12.

On the interval 0. to 1.25000d-01

Label	Value
with weighted error	2.58e-17
log weighted error	16.59
significant figures required	15.22
decimal places required	17.16

Definition at line 176 of file specfunc.cpp.

Referenced by K1(), LPK01(), and LPK1().

◆ ak1_cs

double ak1_cs[17]

static

Initial value:

= {
    2.5 + 0.27443134069738830, 0.07571989953199368,  -0.00144105155647540,
    0.00006650116955125,       -0.00000436998470952, 0.00000035402774997,
    -0.00000003311163779,      0.00000000344597758,  -0.00000000038989323,
    0.00000000004720819,       -0.00000000000604783, 0.00000000000081284,
    -0.00000000000011386,      0.00000000000001654,  -0.00000000000000248,
    0.00000000000000038,       -0.00000000000000006}

Series for ak1.

On the interval 1.25000d-01 to 5.00000d-01

Label	Value
with weighted error	6.06e-17
log weighted error	16.22
significant figures required	15.41
decimal places required	16.83

Definition at line 159 of file specfunc.cpp.

Referenced by K1(), LPK01(), and LPK1().

◆ bi0_cs

double bi0_cs[12]

static

Initial value:

= {
    5.5 - .07660547252839144951, 1.92733795399380827000, .22826445869203013390,
    .01304891466707290428,       .00043442709008164874,  .00000942265768600193,
    .00000014340062895106,       .00000000161384906966,  .00000000001396650044,
    .00000000000009579451,       .00000000000000053339,  .00000000000000000245}

Series for bi0.

On the interval 0. to 9.00000d+00

Label	Value
with weighted error	2.46e-18
log weighted error	17.61
significant figures required	17.90
decimal places required	18.15

Definition at line 126 of file specfunc.cpp.

Referenced by K0(), LPK0(), and LPK01().

◆ bi1_cs

double bi1_cs[11]

static

Initial value:

= {
75 - 0.001971713261099859, 0.407348876675464810, 0.034838994299959456,
001545394556300123,        0.000041888521098377, 0.000000764902676483,
000000010042493924,        0.000000000099322077, 0.000000000000766380,
000000000000004741,        0.000000000000000024}

Series for bi1.

On the interval 0. to 9.00000d+00

Label	Value
with weighted error	2.40e-17
log weighted error	16.62
significant figures required	16.23
decimal places required	17.14

Definition at line 195 of file specfunc.cpp.

Referenced by K1(), LPK01(), and LPK1().

◆ bk0_cs

double bk0_cs[11]

static

Initial value:

= {
    -.5 - 0.03532739323390276872, 0.3442898999246284869,
03597993651536150163,       0.00126461541144692592,
00002286212103119451,       0.00000025347910790261,
00000000190451637722,       0.00000000001034969525,
00000000000004259816,       0.00000000000000013744,
00000000000000000035}

Series for bk0.

On the interval 0. to 4.00000d+00

Label	Value
with weighted error	3.57e-19
log weighted error	18.45
significant figures required	17.99
decimal places required	18.97

Definition at line 73 of file specfunc.cpp.

Referenced by K0(), LPK0(), and LPK01().

◆ bk1_cs

double bk1_cs[11]

static

Initial value:

= {
    1.5 + 0.0253002273389477705, -0.3531559607765448760, -0.1226111808226571480,
    -0.0069757238596398643,      -0.0001730288957513052, -0.0000024334061415659,
    -0.0000000221338763073,      -0.0000000001411488392, -0.0000000000006666901,
    -0.0000000000000024274,      -0.0000000000000000070}

Series for bk1.

On the interval 0. to 4.00000d+00

Label	Value
with weighted error	7.02e-18
log weighted error	17.15
significant figures required	16.73
decimal places required	17.67

Definition at line 144 of file specfunc.cpp.

Referenced by K1(), LPK01(), and LPK1().

◆ hzeta_c

double const hzeta_c[15]

static

Initial value:

= {
    1.00000000000000000000000000000,     0.083333333333333333333333333333,
    -0.00138888888888888888888888888889, 0.000033068783068783068783068783069,
    -8.2671957671957671957671957672e-07, 2.0876756987868098979210090321e-08,
    -5.2841901386874931848476822022e-10, 1.3382536530684678832826980975e-11,
    -3.3896802963225828668301953912e-13, 8.5860620562778445641359054504e-15,
    -2.1748686985580618730415164239e-16, 5.5090028283602295152026526089e-18,
    -1.3954464685812523340707686264e-19, 3.5347070396294674716932299778e-21,
    -8.9535174270375468504026113181e-23}

Coefficients for Maclaurin summation in hzeta().

Evaluated as inverse numbers, i.e. \( \displaystyle\frac{B_{2j}}{(2j)!} \).

Definition at line 205 of file specfunc.cpp.

Referenced by hzeta().

Functions

Variables

Detailed Description

Function Documentation

◆ hzeta()

◆ K0()

◆ K1()

◆ LPK0()

◆ LPK01()

◆ LPK1()

Variable Documentation

◆ ak01_orders

◆ ak02_cs

◆ ak0_cs

◆ ak12_cs

◆ ak1_cs

◆ bi0_cs

◆ bi1_cs

◆ bk0_cs

◆ bk1_cs

◆ hzeta_c