Special Functions

The C code follows the C99 standard, with the usual math functions, as defined in OpenCL. These are automatically available for your models.

Standard functions and constants include:

M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E:

$\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler’s constant $e$

exp, log, pow(x,y), expm1, log1p, sqrt, cbrt:

Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$, $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x) are accurate across all $x$, including $x$ very close to zero.

sin, cos, tan, asin, acos, atan:

Trigonometry functions and inverses, operating on radians.

sinh, cosh, tanh, asinh, acosh, atanh:

Hyperbolic trigonometry functions.

atan2(y,x):

Angle from the $x$-axis to the point $(x,y)$, which is equal to $\tan^{-1}(y/x)$ corrected for quadrant. That is, if $x$ and $y$ are both negative, then atan2(y,x) returns a value in quadrant III where atan(y/x) would return a value in quadrant I. Similarly for quadrants II and IV when $x$ and $y$ have opposite sign.

fabs(x), fmin(x,y), fmax(x,y), trunc, rint:

Floating point functions. rint(x) returns the nearest integer.

NAN:

NaN, Not a Number, $0/0$. Use isnan(x) to test for NaN. Note that you cannot use x == NAN to test for NaN values since that will always return false. NAN does not equal NAN! The alternative, x != x may fail if the compiler optimizes the test away.

INFINITY:

$\infty, 1/0$. Use isinf(x) to test for infinity, or isfinite(x) to test for finite and not NaN.

erf, erfc, tgamma, lgamma: do not use

Special functions that should be part of the standard, but are missing or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma and sas_lgamma instead (see below).

Some non-standard constants and functions are also provided:

M_PI_180, M_4PI_3:

$\frac{\pi}{180}$, $\frac{4\pi}{3}$

SINCOS(x, s, c):

Macro which sets s=sin(x) and c=cos(x). The variables c and s must be declared first.

square(x):

$x^2$

cube(x):

$x^3$

clip(a, a_min, a_max):

$\min(\max(a, a_\text{min}), a_\text{max})$, or NaN if $a$ is NaN.

sas_sinx_x(x):

$\sin(x)/x$, with limit $\sin(0)/0 = 1$.

powr(x, y):

$x^y$ for $x \ge 0$; this is faster than general $x^y$ on some GPUs.

pown(x, n):

$x^n$ for $n$ integer; this is faster than general $x^n$ on some GPUs.

FLOAT_SIZE:

The number of bytes in a floating point value. Even though all variables are declared double, they may be converted to single precision float before running. If your algorithm depends on precision (which is not uncommon for numerical algorithms), use the following:

#if FLOAT_SIZE>4
... code for double precision ...
#else
... code for single precision ...
#endif

SAS_DOUBLE:

A replacement for double so that the declared variable will stay double precision; this should generally not be used since some graphics cards do not support double precision. There is no provision for forcing a constant to stay double precision.

The following special functions and scattering calculations are defined in sasmodels/models/lib. These functions have been tuned to be fast and numerically stable down to $q=0$ even in single precision. In some cases they work around bugs which appear on some platforms but not others, so use them where needed. Add the files listed in source = ["lib/file.c", ...] to your model.py file in the order given, otherwise these functions will not be available.

polevl(x, c, n):

Polynomial evaluation $p(x) = \sum_{i=0}^n c_i x^i$ using Horner’s method so it is faster and more accurate.

$c = \{c_n, c_{n-1}, \ldots, c_0 \}$ is the table of coefficients, sorted from highest to lowest.

source = ["lib/polevl.c", ...] (link to code)

p1evl(x, c, n):

Evaluation of normalized polynomial $p(x) = x^n + \sum_{i=0}^{n-1} c_i x^i$ using Horner’s method so it is faster and more accurate.

$c = \{c_{n-1}, c_{n-2} \ldots, c_0 \}$ is the table of coefficients, sorted from highest to lowest.

source = ["lib/polevl.c", ...] (polevl.c)

sas_gamma(x):

Gamma function sas_gamma$(x) = \Gamma(x)$.

The standard math function, tgamma(x), is unstable for $x < 1$ on some platforms.

source = ["lib/sas_gamma.c", ...] (sas_gamma.c)

sas_gammaln(x):

log gamma function sas_gammaln$(x) = \log \Gamma(|x|)$.

The standard math function, lgamma(x), is incorrect for single precision on some platforms.

source = ["lib/sas_gammainc.c", ...] (sas_gammainc.c)

sas_gammainc(a, x), sas_gammaincc(a, x):

Incomplete gamma function sas_gammainc$(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ and complementary incomplete gamma function sas_gammaincc$(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$

source = ["lib/sas_gammainc.c", ...] (sas_gammainc.c)

sas_erf(x), sas_erfc(x):

Error function sas_erf$(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$ and complementary error function sas_erfc$(x) = \frac{2}{\sqrt\pi}\int_x^{\infty} e^{-t^2}\,dt$.

The standard math functions erf(x) and erfc(x) are slower and broken on some platforms.

source = ["lib/polevl.c", "lib/sas_erf.c", ...] (sas_erf.c)

sas_J0(x):

Bessel function of the first kind sas_J0$(x)=J_0(x)$ where $J_0(x) = \frac{1}{\pi}\int_0^\pi \cos(x\sin(\tau))\,d\tau$.

The standard math function j0(x) is not available on all platforms.

source = ["lib/polevl.c", "lib/sas_J0.c", ...] (sas_J0.c)

sas_J1(x):

Bessel function of the first kind sas_J1$(x)=J_1(x)$ where $J_1(x) = \frac{1}{\pi}\int_0^\pi \cos(\tau - x\sin(\tau))\,d\tau$.

The standard math function j1(x) is not available on all platforms.

source = ["lib/polevl.c", "lib/sas_J1.c", ...] (sas_J1.c)

sas_JN(n, x):

Bessel function of the first kind and integer order $n$, sas_JN$(n, x) =J_n(x)$ where $J_n(x) = \frac{1}{\pi}\int_0^\pi \cos(n\tau - x\sin(\tau))\,d\tau$. If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively.

Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100].

The standard math function jn(n, x) is not available on all platforms.

source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c", ...] (sas_JN.c)

sas_Si(x):

Sine integral Si$(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.

Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100].

This function uses Taylor series for small and large arguments:

For large arguments use the following Taylor series,

$$\text{Si}(x) \sim \frac{\pi}{2} - \frac{\cos(x)}{x}\left(1 - \frac{2!}{x^2} + \frac{4!}{x^4} - \frac{6!}{x^6} \right) - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right)$$

For small arguments,

$$\text{Si}(x) \sim x - \frac{x^3}{3\times 3!} + \frac{x^5}{5 \times 5!} - \frac{x^7}{7 \times 7!} + \frac{x^9}{9\times 9!} - \frac{x^{11}}{11\times 11!}$$

source = ["lib/Si.c", ...] (Si.c)

sas_3j1x_x(x):

Spherical Bessel form sph_j1c$(x) = 3 j_1(x)/x = 3 (\sin(x) - x \cos(x))/x^3$, with a limiting value of 1 at $x=0$, where $j_1(x)$ is the spherical Bessel function of the first kind and first order.

This function uses a Taylor series for small $x$ for numerical accuracy.

source = ["lib/sas_3j1x_x.c", ...] (sas_3j1x_x.c)

sas_2J1x_x(x):

Bessel form sas_J1c$(x) = 2 J_1(x)/x$, with a limiting value of 1 at $x=0$, where $J_1(x)$ is the Bessel function of first kind and first order.

source = ["lib/polevl.c", "lib/sas_J1.c", ...] (sas_J1.c)

Gauss76Z[i], Gauss76Wt[i]:

Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively, computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$.

Similar arrays are available in gauss20.c for 20-point quadrature and in gauss150.c for 150-point quadrature. The macros GAUSS_N, GAUSS_Z and GAUSS_W are defined so that you can change the order of the integration by selecting an different source without touching the C code.

source = ["lib/gauss76.c", ...] (gauss76.c)

Complex numbers

A small complex number library is available using source = ["lib/cl_complex.h", ...]. Numbers are defined as type cdouble, which can be passed to and returned from functions. Operations are as follows:

declare: cdouble x

define: cplx(real, imag)

x.real: creal(x)

x.imag: cimag(x)

1j: I

real + x: radd(real, x)

real - x: rsub(real, x)

real * x: rmul(real, x)

real / x: rdiv(real, x)

x + real: radd(real, x)

x - real: radd(-real, x)

x * real: rmul(real, x)

x / real: rmul(1.0/real, x)

x + y: cadd(x, y)

x - y: csub(x, y)

x * y: cmul(x, y)

x / y: cdiv(x, y)

-x: cneg(x)

abs(x): cabs(x)

angle(x): carg(x)

special functions:

csqrt, cexp, cpow, clog, clog10

csin, ccos, ctan

csinh, ccosh, ctanh

Python Functions

To ease the transition between python and C models, the Special Functions available in C models are reproduced in Python with the same names. All of them are available for import from sasmodels.special. For example, you can include the spherical form factor sas_3j1x_x using from sasmodels.special import sas_3j1x_x in your python model. The python functions operate on both scalars and arrays.

Unlike C models where the constant FLOAT_SIZE is 4 or 8 depending on whether the model is single or double precision, sasmodels.special.FLOAT_SIZE is fixed at 8. Python models should check input.dtype == np.float32 or input.dtype == np.float64 if needed.

SAS_DOUBLE is for C type definitions; it is not available in Python.

sas_j1 = (sin(x) - x cos(x))/x**2 is an additional function that is not available for C.

For Gauss-Lobatto integration points use from sasmodels.special import gauss76 as gauss, then access them with attributes gauss.n, gauss.z and gauss.w. This is equivalent to the macros GAUSS_N, GAUSS_Z and GAUSS_W in C models.

The complex number support functions aren’t needed in Python, and haven’t been implemented.

WARNING: The python functions do not yet have the same accuracy as the corresponding C functions.