Skip to content

Math Operations

Jump to bottom
rizukirr edited this page Jun 11, 2026 · 4 revisions

Math Operations

numc provides element-wise arithmetic, unary math functions, comparisons, reductions, and matrix multiplication. All operations work on all 10 data types and support both contiguous and non-contiguous (transposed, sliced) arrays.

Operation Pattern

Most operations follow the same pattern:

int result = numc_<op>(input_arrays..., output_array);
  • The output array must be pre-allocated with the correct shape and dtype.
  • Returns 0 on success, a negative NumcErrorCode on failure.
  • Input and output arrays must have the same dtype (no implicit promotion).

Element-Wise Binary Operations

Binary operations apply to each corresponding pair of elements. Inputs must have compatible shapes (same shape, or broadcastable -- see Broadcasting).

size_t shape[] = {2, 3};
NumcArray *a = numc_array_create(ctx, shape, 2, NUMC_DTYPE_FLOAT32);
NumcArray *b = numc_array_create(ctx, shape, 2, NUMC_DTYPE_FLOAT32);
NumcArray *out = numc_array_create(ctx, shape, 2, NUMC_DTYPE_FLOAT32);

float da[] = {10, 20, 30, 40, 50, 60};
float db[] = {1, 2, 3, 4, 5, 6};
numc_array_write(a, da);
numc_array_write(b, db);

numc_add(a, b, out);   // [11, 22, 33, 44, 55, 66]
numc_sub(a, b, out);   // [9, 18, 27, 36, 45, 54]
numc_mul(a, b, out);   // [10, 40, 90, 160, 250, 360]
numc_div(a, b, out);   // [10, 10, 10, 10, 10, 10]
Function Formula Notes
numc_add(a, b, out) out = a + b
numc_sub(a, b, out) out = a - b
numc_mul(a, b, out) out = a * b
numc_div(a, b, out) out = a / b Integer division truncates
numc_pow(a, b, out) out = a ^ b
numc_maximum(a, b, out) out[i] = max(a[i], b[i]) Per-element maximum
numc_minimum(a, b, out) out[i] = min(a[i], b[i]) Per-element minimum
numc_fma(a, b, c, out) out = a * b + c Fused multiply-add

Example: Element-Wise Maximum

float da[] = {1, 5, 3, 8, 2, 7};
float db[] = {4, 2, 6, 1, 9, 3};
// ... (create and write arrays)
numc_maximum(a, b, out);   // [4, 5, 6, 8, 9, 7]
numc_minimum(a, b, out);   // [1, 2, 3, 1, 2, 3]

Example: Fused Multiply-Add

// out = a * b + c  (computed in one pass, better precision)
numc_fma(a, b, c, out);

Scalar Operations

Operate on each element with a scalar double value. The scalar is automatically cast to the array's dtype.

numc_add_scalar(a, 100.0, out);   // out = a + 100
numc_sub_scalar(a, 5.0, out);     // out = a - 5
numc_mul_scalar(a, 0.5, out);     // out = a * 0.5
numc_div_scalar(a, 3.0, out);     // out = a / 3

In-Place Scalar Operations

Modify the array directly without needing an output array:

numc_add_scalar_inplace(a, 1000.0);  // a += 1000
numc_sub_scalar_inplace(a, 5.0);     // a -= 5
numc_mul_scalar_inplace(a, 2.0);     // a *= 2
numc_div_scalar_inplace(a, 10.0);    // a /= 10

Unary Operations

Function Formula In-Place Version
numc_neg(a, out) out = -a numc_neg_inplace(a)
numc_abs(a, out) out = |a| numc_abs_inplace(a)
numc_sqrt(a, out) out = sqrt(a) numc_sqrt_inplace(a)
numc_exp(a, out) out = e^a numc_exp_inplace(a)
numc_log(a, out) out = ln(a) numc_log_inplace(a)
numc_clip(a, out, min, max) out = clamp(a, min, max) numc_clip_inplace(a, min, max)

Example: Unary Operations

float da[] = {1.0f, 2.0f, 4.0f, 8.0f};
// ... create a and out ...

numc_log(a, out);   // [0.0, 0.693, 1.386, 2.079]
numc_exp(a, out);   // [2.718, 7.389, 54.598, 2980.958]
numc_sqrt(a, out);  // [1.0, 1.414, 2.0, 2.828]

Example: Clip

Clamps values to a range:

float da[] = {-5.0f, 0.0f, 3.0f, 10.0f};
numc_clip(a, out, 0.0, 5.0);    // [0, 0, 3, 5]
numc_clip_inplace(a, -1.0, 1.0); // [-1, 0, 1, 1]

Integer Types with Math Functions

exp, log, and sqrt on integer types cast through float internally, compute the result, then truncate back:

// int32: exp([0, 1, 10])
// Results: [1, 2, 22026]  (truncated from [1.0, 2.718, 22026.47])

Activation Functions

Function Formula In-Place Version
numc_tanh(a, out) out = tanh(a) numc_tanh_inplace(a)
numc_sigmoid(a, out) out = 1 / (1 + exp(-a)) numc_sigmoid_inplace(a)

Both activations support all 10 dtypes. Float32/float64 use a numerically-stable split formula and ISA-optimized vector kernels (AVX2, AVX-512, NEON, SVE, RVV). Integer dtypes cast through float, compute, then truncate back — useful for testing pipelines but rarely meaningful in practice.

float da[] = {-2.0f, -0.5f, 0.0f, 0.5f, 2.0f};
// ... create a and out (shape {5}, FLOAT32) ...

numc_tanh(a, out);     // [-0.964, -0.462, 0.0, 0.462, 0.964]
numc_sigmoid(a, out);  // [0.119,  0.378, 0.5, 0.622, 0.881]

// In-place variants overwrite a:
numc_tanh_inplace(a);
numc_sigmoid_inplace(a);

Saturation behavior: both functions are bounded — tanh saturates to ±1, sigmoid saturates to [0, 1]. Extremes (e.g. a = 1e30) produce the asymptote without overflow.

Comparison Operations

Comparisons write 1 (true) or 0 (false) into the output. The two inputs must share a dtype, and out must be uint8 — passing any other output dtype returns NUMC_ERR_TYPE. The uint8 result can be fed directly to numc_where as a condition mask (see below).

NumcArray *mask = numc_array_zeros(ctx, shape, ndim, NUMC_DTYPE_UINT8);
numc_gt(a, b, mask);   // mask[i] = (a[i] > b[i])
Function Formula
numc_eq(a, b, out) out[i] = (a[i] == b[i])
numc_gt(a, b, out) out[i] = (a[i] > b[i])
numc_lt(a, b, out) out[i] = (a[i] < b[i])
numc_ge(a, b, out) out[i] = (a[i] >= b[i])
numc_le(a, b, out) out[i] = (a[i] <= b[i])

Scalar Comparisons

numc_eq_scalar(a, 0.0, out);   // out[i] = (a[i] == 0)
numc_gt_scalar(a, 5.0, out);   // out[i] = (a[i] > 5)
numc_lt_scalar(a, 0.0, out);   // out[i] = (a[i] < 0)
numc_ge_scalar(a, 1.0, out);   // out[i] = (a[i] >= 1)
numc_le_scalar(a, 10.0, out);  // out[i] = (a[i] <= 10)

Conditional Selection: numc_where

Selects elements from two arrays based on a condition:

// out[i] = cond[i] ? a[i] : b[i]
numc_where(cond, a, b, out);

a, b, and out must share a dtype. cond must be either that same dtype or uint8 — the uint8 case lets a comparison mask drive the selection directly:

NumcArray *mask = numc_array_zeros(ctx, shape, ndim, NUMC_DTYPE_UINT8);
numc_gt(a, b, mask);          // mask = (a > b), uint8
numc_where(mask, a, b, out);  // out[i] = mask[i] ? a[i] : b[i]

Any nonzero value in cond is treated as true.

Reductions

Reductions collapse one or more dimensions of an array.

Full Reductions

Reduce the entire array to a single value:

size_t shape[] = {2, 3};
NumcArray *a = numc_array_create(ctx, shape, 2, NUMC_DTYPE_FLOAT32);
float da[] = {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
numc_array_write(a, da);

// Output must be a 1-element array
NumcArray *scalar = numc_array_zeros(ctx, (size_t[]){1}, 1, NUMC_DTYPE_FLOAT32);

numc_sum(a, scalar);   // 21.0
numc_mean(a, scalar);  // 3.5
numc_max(a, scalar);   // 6.0
numc_min(a, scalar);   // 1.0

Axis Reductions

Reduce along a specific dimension:

// a is shape (2, 3):
// [[1, 2, 3],
//  [4, 5, 6]]

// Sum along axis 0 (reduce rows): shape (3,)
NumcArray *out0 = numc_array_zeros(ctx, (size_t[]){3}, 1, NUMC_DTYPE_FLOAT32);
numc_sum_axis(a, 0, /*keepdim=*/0, out0);
// Result: [5, 7, 9]

// Sum along axis 1 (reduce columns): shape (2,)
NumcArray *out1 = numc_array_zeros(ctx, (size_t[]){2}, 1, NUMC_DTYPE_FLOAT32);
numc_sum_axis(a, 1, /*keepdim=*/0, out1);
// Result: [6, 15]

keepdim Parameter

When keepdim=1, the reduced axis is kept with size 1 instead of removed:

// Sum axis=0, keepdim=1: shape (2,3) -> (1,3)
NumcArray *kd = numc_array_zeros(ctx, (size_t[]){1, 3}, 2, NUMC_DTYPE_FLOAT32);
numc_sum_axis(a, 0, /*keepdim=*/1, kd);
// Result: [[5, 7, 9]]

This is useful for broadcasting the result back against the original array.

3D Axis Reductions

Works on any number of dimensions:

// b is shape (2, 3, 4)
// Sum along axis 1: shape (2, 3, 4) -> (2, 4)
NumcArray *out3d = numc_array_zeros(ctx, (size_t[]){2, 4}, 2, NUMC_DTYPE_INT32);
numc_sum_axis(b, 1, 0, out3d);

Available Reduction Functions

Full Reduction Axis Reduction Description
numc_sum(a, out) numc_sum_axis(a, axis, keepdim, out) Sum of elements
numc_mean(a, out) numc_mean_axis(a, axis, keepdim, out) Mean of elements
numc_max(a, out) numc_max_axis(a, axis, keepdim, out) Maximum value
numc_min(a, out) numc_min_axis(a, axis, keepdim, out) Minimum value

Argmax / Argmin

Find the index of the maximum or minimum element. Output must be an integer type (INT64 recommended):

// Full argmax/argmin (flattened index)
NumcArray *idx = numc_array_zeros(ctx, (size_t[]){1}, 1, NUMC_DTYPE_INT64);
numc_argmax(a, idx);   // Index of max element in flattened array
numc_argmin(a, idx);   // Index of min element in flattened array

// Axis argmax/argmin
// a is shape (2, 3) -> argmax along axis=1 gives shape (2,)
NumcArray *aidx = numc_array_zeros(ctx, (size_t[]){2}, 1, NUMC_DTYPE_INT64);
numc_argmax_axis(a, 1, 0, aidx);   // [2, 2] (column index of max in each row)

Matrix Multiplication

numc_matmul

High-performance matrix multiply with SIMD-optimized kernels:

// A (2x3) @ B (3x2) = C (2x2)
size_t sa[] = {2, 3}, sb[] = {3, 2}, sc[] = {2, 2};
NumcArray *a = numc_array_create(ctx, sa, 2, NUMC_DTYPE_FLOAT32);
NumcArray *b = numc_array_create(ctx, sb, 2, NUMC_DTYPE_FLOAT32);
NumcArray *c = numc_array_zeros(ctx, sc, 2, NUMC_DTYPE_FLOAT32);

float da[] = {1, 2, 3, 4, 5, 6};
float db[] = {1, 2, 3, 4, 5, 6};
numc_array_write(a, da);
numc_array_write(b, db);

numc_matmul(a, b, c);
// c = [[22, 28], [49, 64]]

numc automatically selects the fastest kernel based on matrix size:

  1. GEMMSUP for small matrices (unpacked SIMD, avoids packing overhead)
  2. Packed GEMM for large matrices (Goto's 5-loop algorithm with cache blocking)
  3. Naive fallback for non-SIMD platforms

numc_dot

Multi-purpose dot product. Behavior depends on input dimensions:

// 1D * 1D: inner product -> scalar
NumcArray *out = numc_array_zeros(ctx, (size_t[]){1}, 1, NUMC_DTYPE_FLOAT32);
numc_dot(vec_a, vec_b, out);

// 2D * 2D: matrix multiplication (same as numc_matmul)
numc_dot(mat_a, mat_b, out);

// ND * 1D: sum-product over last axis of a and b
numc_dot(matrix, vector, out);

// 0D * ND: scalar multiplication
numc_dot(scalar, matrix, out);

numc_matmul_naive

Triple-loop reference implementation. Slower but useful for validation:

numc_matmul_naive(a, b, c);

Clone this wiki locally