Skip to content

Mathematical Functions and Complex Processing

KFR provides element-wise arithmetic, elementary functions, and complex-number operations that use the same vocabulary for individual values and SIMD data. Include <kfr/base.hpp> for the complete expression-capable interface, or <kfr/math.hpp> when only immediate scalar and vector calculations are needed.

#include <kfr/base.hpp>

using namespace kfr;

This article describes elementary functions that take scalar, complex, or vector arguments. The exact domain differs by function: ordinary trigonometric and logarithmic functions operate on real scalar or vector values, while the c-prefixed functions are the complex equivalents. Vector arguments are processed lane by lane, and scalar operands broadcast naturally when mixed with vectors.

Almost every function listed here also has a lazy expression overload. Pass an expression, such as a univector, counter, or arithmetic expression, and it builds a template expression instead of immediately traversing a buffer, so multi-stage pipelines over large data buffers don't need intermediate arrays. See Expressions for expression lifetime and evaluation rules.

Note

Integer inputs to the real transcendental functions are promoted to a floating-point shape. For example, a vec<i32, 4> passed to sin yields a floating-point vector. Mixed arithmetic and binary math functions use the common value type while preserving scalar, vector, and complex shape where applicable.

Immediate values, vectors, and expressions

The following three forms use the same operations at different evaluation scales:

// A scalar result
const f32 gain = exp(-0.5f);

// Four independent SIMD lanes
const f32x4 phases{ 0.f, c_pi<f32, 1, 2>, c_pi<f32>, c_pi<f32, 3, 2> };
const f32x4 waveform = sin(phases);

// A lazy expression, evaluated when it initializes output
const univector<f32> input = truncate(counter(), 1024) * 0.01f;
const univector<f32> output = tanh(input * 2.f);

Use the immediate form for a single value or an explicit vec. Use the expression form for a pipeline over a container. Both calculate the same per-element mathematical result.

Elementary arithmetic

For arithmetic operators, polynomial evaluation, interpolation, rounding, remainders, saturation, bitwise operations, and numeric conversion, see Arithmetic, Bitwise Operations, and Type Conversion.

Trigonometric functions

sin, cos, and tan take angles in radians. Their degree-based counterparts are sindeg, cosdeg, and tandeg; sindeg(90.f) is the same value as sin(c_pi<f32, 1, 2>).

Inverse trigonometric functions and phase

asin and acos accept real inputs in \([-1,1]\). Their outputs are in \([-\pi/2,\pi/2]\) and \([0,\pi]\), respectively. Inputs outside that domain produce invalid floating-point results.

atan returns an angle in \([-\pi/2,\pi/2]\); atandeg returns the same angle in degrees. atan2 and atan2deg use the conventional argument order (y, x) and select the quadrant from both signs:

const f32 phase_radians = atan2(1.f, 0.f); // pi / 2
const f32 phase_degrees = atan2deg(-1.f, 0.f); // -90

Use atan2 rather than atan when a two-dimensional direction or a complex phase is required, because division before atan loses quadrant information.

Sinc and paired lane operations

sinc uses the unnormalized definition

\[ \operatorname{sinc}(x) = \begin{cases} \sin(x)/x, & |x| > \epsilon, \\ 1, & |x| \leq \epsilon. \end{cases} \]

Its argument is in radians; the DSP-normalized \(\sin(\pi x)/(\pi x)\) definition uses different scaling. The zero case is defined as the limiting value 1.

sincos and cossin produce interleaved vector results rather than scalar pairs. sincos places sin(x) in even lanes and cos(x) in odd lanes; cossin reverses that order. sincosdeg and cossindeg apply the same rule to degree inputs.

const f32x4 paired_angles{ 0.f, 0.f, c_pi<f32, 1, 2>, c_pi<f32, 1, 2> };
const auto paired_values = sincos(paired_angles);
// { sin(0), cos(0), sin(pi/2), cos(pi/2) } = { 0, 1, 1, 0 }

When sin(x) and cos(x) are already available, sin2x, sin3x, cos2x, and cos3x derive multiple-angle values without evaluating trigonometric functions again.

Fast trigonometric approximations

fastsin, fastcos, fastsindeg, and fastcosdeg trade accuracy for speed. Use the regular functions unless reduced precision is appropriate for the algorithm. Measured over \([0, 2\pi]\), the errors are approximately the same for float and double:

Function Maximum absolute error Mean absolute error
fastsin \(1.216 \times 10^{-6}\) \(4.738 \times 10^{-7}\)
fastcos \(1.222 \times 10^{-6}\) \(4.550 \times 10^{-7}\)

Outside \([0, 2\pi]\), these approximations are not numerically stable and their error increases rapidly.

Exponential, logarithmic, and root functions

exp, exp2, and exp10 calculate \(e^x\), \(2^x\), and \(10^x\). exp maps negative infinity to zero and positive infinity to positive infinity. exp_fmadd evaluates \(\exp(xm+a)\) directly.

log, log2, and log10 calculate real logarithms. Use a strictly positive argument: zero produces negative infinity, and a negative argument produces NaN. logb returns the unbiased binary exponent and returns negative infinity for zero. logn calculates a logarithm with a specified base, logm multiplies a natural logarithm by a scale, and log_fmadd evaluates \(\log(x)m+a\).

const f32x4 amplitude{ 1.f, 0.5f, 0.f, -0.5f };
const auto safe_amplitude = max(amplitude, 1.0e-20f);
const auto level = log10(safe_amplitude) * 20.f;

pow calculates a real power. A negative base is valid only with an integer exponent; a negative base combined with a non-integer exponent yields NaN. root calculates a real \(n\)th root. It accepts negative values for odd integer degrees and returns a negative root; even or non-integer roots of negative values yield NaN. sqrt is the positive square root and produces NaN for negative real inputs. cbrt supports negative inputs.

const f32 sixth_power = pow(2.f, 6.f); // 64
const f32 negative_root = root(-32.f, 5.f); // -2
const f32 cube_root = cbrt(-8.f); // -2

Hyperbolic and special functions

The hyperbolic family comprises sinh, cosh, tanh, and coth. coth is undefined at zero. The interleaved vector functions sinhcosh and coshsinh place the first named result in even lanes and the second in odd lanes, following the same convention as sincos.

gamma is an approximation of the gamma function for positive real arguments. factorial_approx evaluates an approximate factorial as \(\Gamma(x+1)\) and is useful for floating-point calculation. factorial is the integer lookup form, returning a uint64_t result for its supported range, including factorial(0) == 1, and a maximum uint64_t sentinel outside it. Values larger than \(20!\) do not fit in uint64_t.

modzerobessel calculates the modified zeroth-order Bessel function of the first kind, \(I_0(x)\). It is commonly useful in coefficient calculations such as Kaiser windows.

Note

gamma, factorial_approx, and modzerobessel are approximate implementations. Exact ulp error bounds and special-value behavior depend on the target version and input range.

Complex numbers

complex is KFR's complex value type, normally an alias of std::complex. c32 and c64 are convenient single- and double-precision aliases. Construct a scalar complex value or a vector of complex values with make_complex. real and imag extract its components, and cconj returns its conjugate.

const c32 z = make_complex(3.f, 4.f);
const f32 re = kfr::real(z); // 3
const f32 im = kfr::imag(z); // 4
const c32 conjugate = cconj(z); // 3 - 4i

For a vec<complex<f32>, N>, each vector lane is one complex number. Its flattened scalar layout is interleaved as re0, im0, re1, im1, .... cdecom exposes that flattened real vector and ccomp rebuilds the complex vector. cdupreal, cdupimag, cswapreim, cnegreal, and cnegimag are SIMD layout transforms for specialized complex algorithms.

Complex magnitude, phase, and coordinates

cabssqr returns \(\operatorname{re}(z)^2 + \operatorname{im}(z)^2\); cabs returns its square root. carg returns the phase using atan2(imag(z), real(z)). polar represents a complex value as another complex value whose real part is magnitude and imaginary part is phase. cartesian reverses that conversion.

const c32 z = make_complex(3.f, 4.f);
const f32 magnitude = cabs(z); // 5
const f32 phase = carg(z);
const c32 magnitude_phase = polar(z);
const c32 restored = cartesian(magnitude_phase);

Complex elementary functions

Use the complex-prefixed functions when the input or desired result is complex:

Operation Functions
Trigonometric and hyperbolic csin, ccos, csinh, ccosh
Exponential cexp, cexp2, cexp10
Logarithm clog, clog2, clog10
Roots and powers csqrt, csqr
Magnitude and phase cabssqr, cabs, carg, polar, cartesian

For example, the complex exponential maps a logarithmic magnitude and phase back to Cartesian form:

const c32 exponent = make_complex(log(2.f), c_pi<f32, 1, 2>);
const c32 value = cexp(exponent); // approximately 0 + 2i

clog uses magnitude and principal phase, so complex logarithms and roots follow their standard branch-dependent interpretation. If a signal-processing algorithm needs phase continuity across a sequence, unwrap or otherwise manage phase at the algorithm level rather than assuming a scalar complex function will do so.

Choosing the right operation

  • Use pow2 through pow5, sqr, or cub when the exponent is fixed and small; use ipow for a non-negative integer exponent and pow for a real exponent.
  • Use atan2 for directions and phases, and remember its (y, x) order.
  • Protect real logarithm inputs with max or another domain policy when input data can reach zero or become negative.
  • Use cabs and carg for complex magnitude and phase rather than manually dividing components.
  • Prefer the lazy expression overloads for multi-stage processing of large buffers; materialize only at an intended output or evaluation boundary.

See Also