special Namespace Reference

Name space for special functions. More...

Functions
template<typename T >
constexpr complex_type< T >	i (static_cast< T >(0), static_cast< T >(1))
	I constant.
template<typename T >
complex_type< T >	erfcx (complex_type< T > z)
	Compute erfcx(z) = exp(z^2)*erfz(z)
template<typename T >
T	sq (T x)
	x^2
template<typename T >
T	w_im_y100 (T x)
	Compute a scaled Dawson integral.
template<typename T >
T	isqpi ()
	Define 1.0/sqrt(pi)
template<typename T >
T	w_im (T x)
	Compute specal case of w(z) = exp(z^2)*erfz(z).
template<typename T >
T	erfcx_y100 (const T y100)
	erfcx(x) = exp(x^2)*erfc(x) function, for real x
template<typename T >
T	erfcx (T x)
	erfcx(x) = exp(x^2)*erfc(x) function, for real x
template<typename T >
T	sinh_taylor (T x)
	sinh(x)
template<typename T >
T	sinc (T x, T sinx)
	sinc(x) = sin(x)/x
template<typename T >
complex_type< T >	w (complex_type< T > z)
	Compute w(z) = exp(z^2)*erfz(z).
template<typename T >
complex_type< T >	taylor (const complex_type< T > z, const T mRe_z2, const T mIm_z2)
	Taylor series.
template<typename T >
complex_type< T >	taylor_erfi (const T x, const T y)
	Taylor erfi series.
template<typename T >
complex_type< T >	erf_complex (const complex_type< T > z)
	Compute the error function erf(z)
template<typename T >
complex_type< T >	erfi (const complex_type< T > z)
	erfi(z) = -i erf(iz)

Variables
template<typename T >
constexpr T	expa2n2 []
	Precomputed table of expa2n2[n-1] = exp(-a2*n*n) for double-precision a2 = 0.26865... in w, below.

Detailed Description

Name space for special functions.

Function Documentation

erf_complex()

template<typename T >

complex_type< T > special::erf_complex

(

const complex_type< T >

z

)

Compute the error function erf(z)

Using w_of_z except for certain regions.

Template Parameters

T	Base type of the calculation.

Parameters

[in]

z

The complex argument.

Returns

erf(z)

erfcx() [1/2]

template<typename T >

complex_type< T > special::erfcx

(

complex_type< T >

z

)

Compute erfcx(z) = exp(z^2)*erfz(z)

Template Parameters

T	Base type of the calculation.

Parameters

[in]

z

The complex argument.

Returns

erfcx(z) = exp(z^2)*erfz(z)

erfcx() [2/2]

template<typename T >

T special::erfcx

(

T

x

)

erfcx(x) = exp(x^2)*erfc(x) function, for real x

Template Parameters

T	Base type of the calculation.

Parameters

[in]

x

The argument.

Returns

erfcx(x) = exp(x^2)*erfc(x)

erfcx_y100()

template<typename T >

T special::erfcx_y100

(

const T

y100

)

erfcx(x) = exp(x^2)*erfc(x) function, for real x

Written by Steven G. Johnson, October 2012.

This function combines a few different ideas.

First, for x > 50, it uses a continued-fraction expansion (same as for the Faddeeva function, but with algebraic simplifications for z=i*x).

Second, for 0 <= x <= 50, it uses Chebyshev polynomial approximations, but with two twists:

a) It maps x to y = 4 / (4+x) in [0,1]. This simple transformation, inspired by a similar transformation in the octave-forge/specfun erfcx by Soren Hauberg, results in much faster Chebyshev convergence than other simple transformations I have examined.

b) Instead of using a single Chebyshev polynomial for the entire [0,1] y interval, we break the interval up into 100 equal subintervals, with a switch/lookup table, and use much lower degree Chebyshev polynomials in each subinterval. This greatly improves performance in my tests.

For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x), with the usual checks for overflow etcetera.

Performance-wise, it seems to be substantially faster than either the SLATEC DERFC function [or an erfcx function derived therefrom] or Cody's CALERF function (from netlib.org/specfun), while retaining near machine precision in accuracy. *‍/

Given y100=100*y, where y = 4/(4+x) for x >= 0, compute erfc(x).

Uses a look-up table of 100 different Chebyshev polynomials for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated with the help of Maple and a little shell script. This allows the Chebyshev polynomials to be of significantly lower degree (about 1/4) compared to fitting the whole [0,1] interval with a single polynomial.

Template Parameters

T	Base type of the calculation.

Parameters

[in]

y100

Interval argument.

Returns

erfcx(x) = exp(x^2)*erfc(x) function, for real x

erfi()

template<typename T >

complex_type< T > special::erfi

(

const complex_type< T >

z

)

erfi(z) = -i erf(iz)

Template Parameters

T	Base type of the calculation.

Parameters

[in]

z

Complex argument.

Returns

erfi(z)

isqpi()

template<typename T >

T special::isqpi

(

)

Define 1.0/sqrt(pi)

Returns

1/sqrt(pi)

sinc()

template<typename T >

T special::sinc	(	T	x,
		T	sinx
	)

sinc(x) = sin(x)/x

Return sinc(x) = sin(x)/x, given both x and sin(x) Note: Since we only use this in cases where sin(x) has already been computed.

Template Parameters

T	Base type of the calculation.

Parameters

[in]	x	Argument of sine.
[in]	sinx	Precomputed sine.

Returns

sinc(x) = sin(x)/x

sinh_taylor()

template<typename T >

T special::sinh_taylor

(

T

x

)

sinh(x)

sinh(x) via Taylor series, accurate to machine precision for |x| < 1E-2

Template Parameters

T	Base type of the calculation.

Parameters

[in]

x

Argument.

Returns

sinh(x)

sq()

template<typename T >

T special::sq

(

T

x

)

x^2

Template Parameters

T	Base type of the calculation.

Parameters

[in]

x

Argument.

Returns

x^2

taylor()

template<typename T >

complex_type< T > special::taylor	(	const complex_type< T >	z,
		const T	mRe_z2,
		const T	mIm_z2
	)

Taylor series.

Use Taylor series for small |z|, to avoid cancellation inaccuracy erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)

Template Parameters

T	Base type of the calculation.

Parameters

[in]	z	Complex argument.
[in]	mRe_z2	Real argument squared.
[in]	mIm_z2	Imaginary argument squared.

Returns

erf(z) where z = x + Iy.

taylor_erfi()

template<typename T >

complex_type< T > special::taylor_erfi	(	const T	x,
		const T	y
	)

Taylor erfi series.

For small |x| and small |xy|, use Taylor series to avoid cancellation inaccuracy:

erf(x+iy) = erf(iy) + 2*exp(y^2)/sqrt(pi) * [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ...

• i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]

where: erf(iy) = exp(y^2)Im[w(y)]

Template Parameters

T	Base type of the calculation.

Parameters

[in]	x	Real argument.
[in]	y	imaginary argument.

Returns

erf(x+iy)

w()

template<typename T >

complex_type< T > special::w

(

complex_type< T >

z

)

Compute w(z) = exp(z^2)*erfz(z).

Template Parameters

T	Base type of the calculation.

Parameters

[in]

z

The complex argument.

Returns

erfcx(z) = exp(z^2)*erfz(z)

w_im()

template<typename T >

T special::w_im

(

T

x

)

Compute specal case of w(z) = exp(z^2)*erfz(z).

Specialized for imaginary arguments.

Template Parameters

T	Base type of the calculation.

Parameters

[in]

x

The argument.

Returns

erfcx(z) = exp(z^2)*erfz(z)

w_im_y100()

template<typename T >

T special::w_im_y100

(

T

x

)

Compute a scaled Dawson integral.

w_im(x) = 2*Dawson(x)/sqrt(pi) equivalent to the imaginary part w(x) for real x.

Uses methods similar to the erfcx calculation above: continued fractions for large |x|, a lookup table of Chebyshev polynomials for smaller |x|, and finally a Taylor expansion for |x|<0.01. Given y100 = 100*y, where y = 1/(1 + x) for x >= 0, compute w_im(x).

Template Parameters

T	Base type of the calculation.

Parameters

[in]

x

The real argument.