# Mathematical Operations

Keywords: math transcendental functions

The math module implements several useful functions for digital signal processing including transcendental function not necessarily in the standard C library, windowing functions, and polynomial manipulation methods.

## Transcendental Functions∞

This section describes the implementation and interface to transcendental functions not in the C standard library including a full arrangement of Gamma and Bessel functions.[ref:tab-math-transcendentals] summarizes the interfaces provided in liquid .

Table [tab-math-transcendentals]. Summary of Transcendental Math Interfaces

 Function Interface $$\ln\Gamma(z)$$ liquid_lngammaf(z) $$\Gamma(z)$$ liquid_gammaf(z) $$\ln\gamma(z,\alpha)$$ liquid_lnlowergammaf(z,alpha) $$\gamma(z,\alpha)$$ liquid_lowergammaf(z,alpha) $$\ln\Gamma(z,\alpha)$$ liquid_lnuppergammaf(z,alpha) $$\Gamma(z,\alpha)$$ liquid_uppergammaf(z,alpha) $$n!$$ liquid_factorialf(n) $$\ln I_\nu(z)$$ liquid_lnbesselif(nu,z) $$I_\nu(z)$$ liquid_besselif(nu,z) $$I_0 (z)$$ liquid_besseli0f(z) $$J_\nu(z)$$ liquid_besseljf(nu,z) $$J_0 (z)$$ liquid_besselj0f(z) $$Q(z)$$ liquid_Qf(z) $$Q_M(\alpha,\beta)$$ liquid_MarcumQf(M,alpha,beta) $$Q_1(\alpha,\beta)$$ liquid_MarcumQ1f(alpha,beta) $$\sinc(z)$$ liquid_sincf(z) $$\lceil\log_2(n)\rceil$$ liquid_nextpow2(n) $${n \choose k}$$ liquid_nchoosek(n,k)

### liquid_gammaf(z), liquid_lngammaf(z)∞

liquid computes $$\Gamma(z)$$ from $$\ln\Gamma(z)$$ (see below) due to its steep, exponential response to $$z$$ . The complete Gamma function is defined as

$$\Gamma(z) \triangleq \int_0^\infty{t^{z-1}e^{-t}dt}$$

The upper an lower incomplete Gamma functions are described in Sections[ref:section-math-transcendentals-uppergamma] and[ref:section-math-transcendentals-lowergamma] , respectively. The natural log of the complete Gamma function is computed by splitting into discrete piecewise sections:

$$\ln\left[ \Gamma(z) \right] \approx \begin{cases} % undefined for z \lt 0 \text{undefined} & z \lt 0 \\ % low signal approximation (recursion) \ln\Gamma(z+1) - \ln(z) & 0 \le z \lt 10 \\ % high signal approximation \frac{z}{2} \ln\left( \frac{2\pi}{z} \right) \left( \ln\left(z + \frac{1}{12 z - 0.1/z} \right) - 1 \right) & z \ge 0.6 \end{cases}$$

### liquid_lowergammaf(z,a), liquid_lnlowergammaf(z,a) (lower incomplete Gamma)∞

Like $$\Gamma(z)$$ , liquid computes the lower incomplete gamma function$$\gamma(z,\alpha)$$ from its logarithm $$\ln\gamma(z,\alpha)$$ due to its steep, exponential response to $$z$$ . The lower incomplete Gamma function is defined as

$$\gamma(z,\alpha) \triangleq \int_0^\alpha{ t^{z-1}e^{-t}dt }$$

liquid computes the log of lower incomplete Gamma function as

$$\ln\gamma(z,\alpha) = z \ln(\alpha) + \ln\Gamma(z) - \alpha + \ln\left[ \sum_{k=0}^\infty{ \frac{\alpha^k}{\Gamma(z + k + 1)} } \right]$$

### liquid_uppergammaf(z,a), liquid_lnuppergammaf(z,a) (upper incomplete Gamma)∞

Like $$\Gamma(z)$$ , liquid computes the upper incomplete gamma function$$\Gamma(z,\alpha)$$ from $$\ln\Gamma(z,\alpha)$$ due to its steep, exponential response to $$z$$ . The complete Gamma function is defined as

$$\Gamma(z,\alpha) \triangleq \int_\alpha^\infty{t^{z-1}e^{-t}dt}$$

By definition the sum of the lower and upper incomplete gamma functions is the complete Gamma function:$$\Gamma(z) = \gamma(z,\alpha) + \Gamma(z,\alpha)$$ . As such, liquid computes the upper incomplete Gamma function as

$$\Gamma(z,\alpha) = \Gamma(z) - \gamma(z,\alpha)$$

### liquid_factorialf(n)∞

liquid computes $$n!=n\cdot(n-1)\cdot(n-2)\cdots3\cdot2\cdot1$$ iteratively for small values of $$n$$ , and with the Gamma function for larger values. Specifically, $$n! = \Gamma(n+1)$$ .

### liquid_nchoosek()∞

liquid computes binomial coefficients using the liquid_nchoosek() method:

$${n \choose k} = \frac{n!}{(n-k)!k!}$$

Because the arguments can explode for relatively large values of $$n$$ and$$k$$ , liquid uses the following approximation under certain conditions:

$${n \choose k} \approx \exp\Bigl\{ \ln\Gamma(n+1) - \ln\Gamma(n-k+1) - \ln\Gamma(k+1) \Bigr\}$$

### liquid_nextpow2()∞

computes $$\lceil \log_2(x) \rceil$$

### liquid_sinc(z)∞

The $$\sinc$$ function is defined as

$$\sinc(z) = \frac{\sin(\pi z)}{\pi z}$$

Simply evaluating the above equation with finite precision for $$z$$ results in a discontinuity for small $$z$$ , and is approximated by expanding the first few terms of the series

$$\sinc(z) = \prod_{k=1}^{\infty}{ \cos\left( 2^{-k} \pi z \right) }$$

### liquid_lnbesselif(), liquid_besselif(), liquid_besseli0f()∞

$$I_\nu(z)$$ is the modified Bessel function of the first kind and is particularly useful for filter design. An iterative method for computing $$I_\nu$$ comes from Gross(1995),

$$I_\nu(z) = \left(\frac{z}{2}\right)^\nu \sum_{k=0}^{\infty}{\frac{\left(\frac{1}{4}z^2\right)^k}{k!\Gamma(k+\nu+1)}}$$

Due to its steep response to $$z$$ it is often useful to compute$$I_\nu(z)$$ by first computing $$\ln I_\nu(z)$$ as

..eqnarray [eqn-math-lnbesseli]

\ln I_\nu(z) & = &
\nu\ln(z/2) + \ln\left[
\sum_{k=0}^\infty{
\frac{
\left(\frac{1}{4}z^2\right)^k
}{
k! \Gamma(\nu + k + 1)
}
}
\right] \\
& = &
\nu\ln(z/2) + \ln\left[
\sum_{k=0}^\infty{ \exp\Bigl\{
2 k \ln(z/2) - \ln\Gamma(k+1) - \ln\Gamma(\nu+k+1)
\Bigr\} }
\right] \\


For $$\nu=0$$ a good approximation can be derived by using piecewise polynomials,

$$\ln\Bigl[\ln\left(I_0(z)\right)\Bigr] \approx c_0 + c_1 t + c_2 t ^2 + c_3 t^3$$

where $$t=\ln(z)$$ and

$$\left\{c_0,c_1,c_2,c_3\right\} = \begin{cases} \left\{\text{-1.52624, 1.9597, -9.4287e-03, -7.0471e-04}\right\} & t \lt 0.5 \\ \left\{\text{-1.5531, 1.8936, -0.07972, -0.01333}\right\} & 0.5 \le t \lt 2.3 \\ \left\{\text{-1.2958, 1.7693, -0.1175, 0.006341}\right\} & \text{else}. \end{cases}$$

This is a particularly useful approximation for the Kaiser window in fixed-point math where $$w[n]$$ is computed as the ratio of two large numbers.

### liquid_lnbesseljf(), liquid_besselj0f()∞

$$J_\nu(z)$$ is the Bessel function of the first kind and is found in Doppler filter design. liquid computes $$J_\nu(z)$$ using the series expansion

$$J_\nu(z) = \sum_{k=0}^\infty{ \frac{ (-1)^k }{ 2^{2k+|v|} k! \left(|v|+k\right)! } z^{2k+|v|} }$$

### liquid_Qf(), liquid_MarcumQf(), liquid_MarcumQ1f()∞

The $$Q$$ -function is commonly used in signal processing and is defined as

.. eqnarray
Q(z) &=& \frac{1}{2}\left(1 - \textup{erf}(z/\sqrt{2})\right) \\
&=& \frac{1}{\sqrt{2\pi}} \int_{z}^{\infty} { \exp\left\{-u^2/2 \right\} du }


Similarly Marcum's $$Q$$ -function is defined as the following, with an appropriate expansion:

.. eqnarray

Q_M(\alpha,\beta) & = &
\int_{\beta}^{\infty}{
u\left(\frac{u}{\alpha}\right)^{M-1}
\exp\left\{ -\frac{u^2 + \alpha^2}{2}\right\}
I_{M-1}(\alpha u)
du
}\\
& = &
\exp\left\{-\frac{\alpha^2 + \beta^2}{2}\right\}
\sum_{k=1-M}^{\infty}{
\left(\frac{\alpha}{\beta}\right)^k
I_k(\alpha\beta)
}


where $$I_\nu$$ is the modified Bessel function of the first kind (see [ref:section-math-transcendentals-besseli] ). liquid implements $$Q_M(\alpha,\beta)$$ with the function liquid_MarcumQf(M,a,b) using the approximation {cite:Helstrom:1992((25))}

$$Q_M(\alpha,\beta) \approx \textup{erfc}(u), \,\,\, u=\frac{\beta-\alpha-M}{\sigma^2}, \,\,\, \sigma = M + 2\alpha$$

which works over a reasonable range of $$M$$ , $$\alpha$$ , and $$\beta$$ . The special case for $$M=1$$ is implemented in liquid using the function liquid_MarcumQ1f(M,a,b) using the expansion {cite:Helstrom:1960},

$$Q_1(\alpha,\beta) = \exp\left\{-\frac{\alpha^2 + \beta^2}{2}\right\} \sum_{k=0}^{\infty}{ \left(\frac{\alpha}{\beta}\right)^k I_k(\alpha\beta) }$$

which converges quickly with a few iterations.

## Complex Trigonometry∞

This section describes the implementation and interface to complex trigonometric functions not in the C standard library.[ref:tab-math-transcendentals] summarizes the interfaces provided in liquid .

Table [tab-math-complex]. Summary of Complex Trigonometric Math Interfaces

 Function Interface $$\sqrt{z}$$ liquid_csqrtf(z) $$e^{z}$$ liquid_cexpf(z) $$\ln(z)$$ liquid_clogf(z) $$\sin^{-1}(z)$$ liquid_casinf(z) $$\cos^{-1}(z)$$ liquid_cacosf(z) $$\tan^{-1}(z)$$ liquid_catanf(z)

### liquid_csqrtf()∞

The function liquid_csqrtf(z) computes the complex square root of a number

$$\sqrt{z} = \sqrt{\frac{r+a}{2}} + j\text{sgn}\bigl(\Im\{z\}\bigr) \sqrt{\frac{r-a}{2}}$$

where $$r=|z|$$ , $$a=\Re\{z\}$$ , and $$\text{sgn}(t)=t/|t|$$ .

### liquid_cexpf()∞

The function liquid_cexpf(z) computes the complex exponential of a number

$$e^{z} = \exp\bigl\{a\bigr\} \bigl( \cos(b) + j\sin(b) \bigr)$$

where $$a=\Re\{z\}$$ and $$b=\Im\{z\}$$ .

### liquid_clogf()∞

The function liquid_clogf(z) computes the complex natural logarithm of a number.

$$\log(z) = \log(|z|) + j\arg(z)$$

### liquid_cacosf()∞

The function liquid_cacosf(z) computes the complex $$\arccos$$ of a number

$$\arccos(z) = \begin{cases} -j \log\bigl( z + \sqrt{z^2 - 1} \bigr) & \text{sgn}\bigl(\Re\{z\}\bigr) = \text{sgn}\bigl(\Im\{z\}\bigr) \\ -j \log\bigl( z - \sqrt{z^2 - 1} \bigr) & \text{otherwise} \end{cases}$$

### liquid_casinf()∞

The function liquid_casinf(z) computes the complex $$\arcsin$$ of a number

$$\arcsin(z) = \frac{\pi}{2} - \arccos(z)$$

### liquid_catanf()∞

The function liquid_catanf(z) computes the complex $$\arctan$$ of a number

$$\arctan(z) = \frac{j}{2} \log\left( \frac{1-jz}{1+jz} \right)$$

## Modular Arithmetic∞

### liquid_is_prime(n)∞

Returns 1 if $$n$$ is prime, 0 otherwise.

### liquid_factor(n,*factors,*num_factors)∞

Computes all the prime factors of a number $$n$$ in ascending order. Example:

unsigned int factors[LIQUID_MAX_FACTORS];
unsigned int num_factors;
liquid_factor(280, factors, &num_factors);


>> num_factors = 5>> factors = { 2, 2, 2, 5, 7 }

### liquid_unique_factor(n,*factors,*num_factors)∞

Computes all the unique prime factors of a number $$n$$ in ascending order. Example:

unsigned int factors[LIQUID_MAX_FACTORS];
unsigned int num_factors;
liquid_unique_factor(280, factors, &num_factors);


>> num_factors = 3>> factors = { 2, 5, 7 }

### liquid_modpow(base,exp,n)∞

Computes $$c = b^x \pmod n$$

### liquid_primitive_root(n)∞

Finds and returns smallest primitive root of a number $$n$$ .

### liquid_primitive_root_prime(n)∞

Finds and returns smallest primitive root of a prime number $$n$$ .

### liquid_totient(n)∞

Euler's totient function $$\varphi(n)$$ computes the number of positive integers less than or equal to $$n$$ that are relatively prime to $$n$$ . The totient function is computed as

$$\varphi(n) = n \prod_{p|n} {\left(1 - \frac{1}{p}\right)}$$

where the product is computed over the unique factors of $$n$$ . For example, if $$n=24 = 2^3 \cdot 3$$ , then$$\varphi(24) = 24\left(1-\frac{1}{2}\right)\left(1-\frac{1}{3}\right)=8$$ .