Source
stdlib/fixarray/fixarray.hpp
1
// Copyright (c) 2026 BigBrain LLC. MIT-licensed (see LICENSE).2
// Original work; see ACKNOWLEDGMENTS.md for the open-source ideas we build upon.3
// cheatah-deps: ndarray4
#pragma once6
/**7
* @file fixarray.hpp8
* @brief cheatah `fixarray` — fixed-extent arrays (@ref cheatah::fixarray::Fixed): exactly like an9
* @ref cheatah::ndarray::NDArray, only faster.10
*11
* An @ref cheatah::ndarray::NDArray carries its shape at runtime and its elements on the heap, which12
* is what makes it general. When the shape is known at compile time and tiny — a 3-D direction, a13
* 4×4 transform — that generality is the whole cost: a heap allocation, a stride computation and an14
* indirection per operation, to move sixteen floats.15
*16
* @ref cheatah::fixarray::Fixed is the same idea with the shape moved into the type. The extents are17
* template parameters, the elements live inline (a `std::array`, so the value is trivially copyable18
* and sits on the stack or straight inside another struct), the loops have compile-time trip counts19
* and auto-vectorize, and nothing allocates. **These are the types to reach for in high-performance20
* applications — a renderer's transforms, a physics solver's contact frames, a filter's small21
* state** — where the same matrix is built and consumed millions of times a second.22
*23
* Everything else is deliberately the same as `NDArray`: element types are the same @ref24
* cheatah::ndarray::Field, the mathematical index is `(row, column)`, the vocabulary is numpy's25
* (@ref dot, @ref matmul, @ref transpose, @ref determinant, @ref inverse), and results agree26
* elementwise. Reach for `NDArray` when the shape is data; reach for `Fixed` when the shape is a27
* fact about the program.28
*29
* **One deliberate difference: a matrix is stored COLUMN-MAJOR**, where `NDArray` is row-major. The30
* indexing you write is unchanged — `m(row, col)` means what it says, and the constructor still31
* takes elements in reading order — but @ref Fixed::data() hands back columns, not rows. Two reasons,32
* both measured: `m * v` becomes a sum of scaled columns, which is contiguous, vertical, and33
* vectorizes, instead of four horizontal dot products that cost a shuffle network; and the buffer is34
* already in the order graphics APIs (GLSL, SPIR-V, Metal) and GLM expect, so uploading a transform35
* is a copy rather than a transpose. Only reach for `data()` when you mean the raw buffer.36
*37
* ```38
* using namespace cheatah::fixarray;39
* vec3f up{0.0F, 1.0F, 0.0F}; // a 3-vector, 12 bytes, no allocation40
* mat4f m = mat4f::identity(); // a 4x4, 64 bytes — exactly a push constant41
* vec3f v = normalize(cross(up, w)); // numpy's vocabulary, glm's speed42
* ```43
*44
* Rank 1 (a vector) and rank 2 (a matrix) are supported; higher ranks are a mechanical extension of45
* the same storage and are added when a caller needs one.46
*47
* This module is templates only — header-only, nothing is compiled into a library — so the caller's48
* optimization flags apply. The `linalg` module remains the home of the heavy, shape-generic numerics49
* on `NDArray` (LU, QR, SVD, eigen); `Fixed` owns the small closed forms where a general factorization50
* would cost more than the answer.51
*52
* **Performance.** Benchmarked against [GLM](https://github.com/g-truc/glm) over the complete overlap53
* of the two APIs — 160 pairs, every operation, sizes 2/3/4, `float` and `double`, with the outputs54
* verified identical before either is timed — `Fixed` is **faster than or at parity with GLM on every55
* one** (20 faster, 140 at parity, none slower; medians over 9 interleaved repetitions, a win counting56
* only above both 1.15x and 0.25 ns). It wins where structure pays: `mat4f::identity()`57
* 2.69×, `mat4f * mat4f` 1.70×, `mat4f + mat4f` 2.03×, `inverse(mat4d)` 1.38×. No intrinsics — the code is58
* shaped so the compiler vectorizes it. A regression gate (`scripts/bench_gate.sh`) keeps it true.59
* See the @ref performance "Small fixed-size math vs GLM" section for the how and the numbers.60
*/62
#include <array>63
#include <cmath>64
#include <concepts>65
#include <cstddef>66
#include <ostream>67
#include <algorithm>68
#include <iterator>69
#include <limits>70
#include <stdexcept>71
#include <string>72
#include <utility>74
#include "ndarray.hpp"76
namespace cheatah::fixarray {78
/// The product of a pack of extents — a fixed array's element count, `1` for an empty pack.79
/// @tparam Dims the extents.80
template <std::size_t... Dims>81
inline constexpr std::size_t extent_product = (std::size_t{1} * ... * Dims);83
/**84
* Constrains a @ref Fixed to a supported rank: 1 (a vector) or 2 (a matrix). Higher ranks are a85
* mechanical extension of the same storage, added when a caller needs one — the concept is what86
* turns "not yet" into a readable compile error instead of a template-instantiation wall.87
* @tparam Rank the number of extents.88
*/89
template <std::size_t Rank>90
concept SupportedRank = (Rank == 1 || Rank == 2);92
/**93
* A fixed-extent, inline-stored array — an @ref cheatah::ndarray::NDArray whose shape lives in the94
* type. Trivially copyable and allocation-free; a matrix is stored column-major (see the file doc).95
*96
* @tparam T the element type; any @ref cheatah::ndarray::Field, exactly as `NDArray` accepts.97
* @tparam Dims the extents. One extent is a vector, two are a matrix (rows, then columns).98
*/99
template <ndarray::Field T, std::size_t... Dims>100
requires SupportedRank<sizeof...(Dims)> && (((Dims > 0) && ...))101
class Fixed {102
public:103
/// The element type.104
using value_type = T;106
/// The number of extents: 1 for a vector, 2 for a matrix.107
static constexpr std::size_t rank = sizeof...(Dims);108
/// The total number of elements.109
static constexpr std::size_t size = extent_product<Dims...>;110
/// The extents, in order (rows, then columns for a matrix).111
static constexpr std::array<std::size_t, rank> shape{Dims...};113
/// Rows — the first extent (a vector has one row).114
static constexpr std::size_t rows = shape[0];115
/// Columns — the second extent, or 1 for a vector.116
static constexpr std::size_t cols = rank == 2 ? shape[1] : 1;118
/// Every element zero — the additive identity, and what a default-constructed value holds.119
/// @complexity O(size).120
/// @alloc none.121
/// @test Fixarray.DefaultIsZero122
constexpr Fixed() = default;124
/**125
* Construct from exactly @ref size elements, written in READING order: a matrix is given row by126
* row, the way it appears on paper, regardless of how it is stored. Arguments are converted to127
* @p T, so a `vec3f` accepts the doubles a cheatah program computes with.128
* @tparam Args the argument types; each must be convertible to @p T.129
* @param args the elements in reading order; exactly @ref size of them.130
* @complexity O(size).131
* @alloc none.132
* @test Fixarray.MatrixIndexing133
* @crtest FixarrayCompileRun.ConstructAndDot134
*/135
template <class... Args>136
requires(sizeof...(Args) == size) && (std::convertible_to<Args, T> && ...)137
explicit constexpr Fixed(Args... args) {138
const std::array<T, size> reading_order{static_cast<T>(args)...};139
if constexpr (rank == 1) {140
data_ = reading_order;141
} else {142
for (std::size_t r = 0; r < rows; ++r) {143
for (std::size_t c = 0; c < cols; ++c) { data_[c * rows + r] = reading_order[r * cols + c]; }144
}145
}146
}148
/**149
* The square identity: ones on the diagonal, zeros elsewhere.150
* @return the identity matrix.151
* @complexity O(size).152
* @alloc none.153
* @test Fixarray.Identity154
*/155
static constexpr Fixed identity()156
requires(rank == 2 && rows == cols)157
{158
// Built element by element in place. Zeroing the buffer and then poking the diagonal would159
// store every byte twice; this stores each once, and the compiler folds it to a constant.160
// In a square column-major buffer the diagonal is exactly the indices divisible by rows + 1.161
return identity_impl(std::make_index_sequence<size>{});162
}164
/**165
* Every element set to @p value — `filled(0)` is the zero value, `filled(1)` a matrix of ones.166
* @param value the element to repeat.167
* @return the filled array.168
* @complexity O(size).169
* @alloc none.170
* @test Fixarray.Filled171
*/172
static constexpr Fixed filled(T value) {173
Fixed result;174
for (std::size_t i = 0; i < size; ++i) { result.data_[i] = value; }175
return result;176
}178
/**179
* Build an array elementwise: each element `i` of the flat, contiguous buffer is `f(i)`. This is180
* the allocation-free, single-pass way to write a component-wise operation — no default zeroing181
* and no separate copy to overwrite, so a call like `abs` or `min` compiles to one vector pass182
* (`minps`/`maxpd`) rather than two. The index `i` runs over the storage order (column-major for183
* a matrix), which is exactly what an elementwise operation wants.184
* @tparam F a callable `T(std::size_t)`.185
* @param f produces element `i` from its flat index.186
* @return the array whose element `i` is `f(i)`.187
* @complexity O(size).188
* @alloc none.189
* @test Fixarray.FromIndices190
*/191
template <class F>192
static constexpr Fixed from_indices(F&& f) {193
return from_indices_impl(std::forward<F>(f), std::make_index_sequence<size>{});194
}196
/**197
* Element @p i of a vector.198
* @param i the index, `0 <= i < size`.199
* @return a reference to the element.200
* @complexity O(1).201
* @alloc none.202
* @test Fixarray.VectorIndexing203
*/204
constexpr T& operator[](std::size_t i)205
requires(rank == 1)206
{207
return data_[i];208
}210
/// The same, indexed by a scoped `enum class` column label (see @ref ndarray::Subscript): the one211
/// place an enum is spent as an index, so `v[Axis::Z]` reads column Z while `Axis` stays strong212
/// everywhere else.213
/// @tparam Ix the enum index type.214
/// @param i the element to address, named by an enumerator.215
/// @return a reference to the element.216
/// @complexity O(1). @alloc none.217
/// @test Fixarray.EnumIndexingOnVectorsAndMatrices218
template <::cheatah::ndarray::Subscript Ix>219
requires(rank == 1 && std::is_enum_v<Ix>)220
constexpr T& operator[](Ix i) {221
return data_[static_cast<std::size_t>(::cheatah::ndarray::subscript_index(i))];222
}224
/**225
* Element @p i of a vector (read-only).226
* @param i the index, `0 <= i < size`.227
* @return a const reference to the element.228
* @complexity O(1).229
* @alloc none.230
* @test Fixarray.VectorIndexing231
*/232
constexpr const T& operator[](std::size_t i) const233
requires(rank == 1)234
{235
return data_[i];236
}238
/// Read-only element by a scoped `enum class` column label (see @ref ndarray::Subscript).239
/// @tparam Ix the enum index type.240
/// @param i the element to address, named by an enumerator.241
/// @return a const reference to the element.242
/// @complexity O(1). @alloc none.243
/// @test Fixarray.EnumIndexingOnVectorsAndMatrices244
template <::cheatah::ndarray::Subscript Ix>245
requires(rank == 1 && std::is_enum_v<Ix>)246
constexpr const T& operator[](Ix i) const {247
return data_[static_cast<std::size_t>(::cheatah::ndarray::subscript_index(i))];248
}250
/**251
* Element (@p row, @p col) of a matrix. The index is mathematical; the storage is column-major.252
* @param row the row, `0 <= row < rows`.253
* @param col the column, `0 <= col < cols`.254
* @return a reference to the element.255
* @complexity O(1).256
* @alloc none.257
* @test Fixarray.MatrixIndexing258
*/259
constexpr T& operator()(std::size_t row, std::size_t col)260
requires(rank == 2)261
{262
return data_[col * rows + row];263
}265
/// The same, with either index a scoped `enum class` label (see @ref ndarray::Subscript) — a named266
/// row or column of a fixed matrix. Mixed integer/enum is allowed; at least one must be an enum, so267
/// the plain `std::size_t` overload still owns the all-integer call.268
/// @tparam R the row index type. @tparam C the column index type; at least one is an enum.269
/// @param row the row to address. @param col the column to address.270
/// @return a reference to the element.271
/// @complexity O(1). @alloc none.272
/// @test Fixarray.EnumIndexingOnVectorsAndMatrices273
template <::cheatah::ndarray::Subscript R, ::cheatah::ndarray::Subscript C>274
requires(rank == 2 && (std::is_enum_v<R> || std::is_enum_v<C>))275
constexpr T& operator()(R row, C col) {276
return (*this)(static_cast<std::size_t>(::cheatah::ndarray::subscript_index(row)),277
static_cast<std::size_t>(::cheatah::ndarray::subscript_index(col)));278
}280
/**281
* Element (@p row, @p col) of a matrix, read-only. Mathematical index; column-major storage.282
* @param row the row, `0 <= row < rows`.283
* @param col the column, `0 <= col < cols`.284
* @return a const reference to the element.285
* @complexity O(1).286
* @alloc none.287
* @test Fixarray.MatrixIndexing288
*/289
constexpr const T& operator()(std::size_t row, std::size_t col) const290
requires(rank == 2)291
{292
return data_[col * rows + row];293
}295
/// Read-only (@p row, @p col) with either index a scoped `enum class` label (see @ref296
/// ndarray::Subscript).297
/// @tparam R the row index type. @tparam C the column index type; at least one is an enum.298
/// @param row the row to address. @param col the column to address.299
/// @return a const reference to the element.300
/// @complexity O(1). @alloc none.301
/// @test Fixarray.EnumIndexingOnVectorsAndMatrices302
template <::cheatah::ndarray::Subscript R, ::cheatah::ndarray::Subscript C>303
requires(rank == 2 && (std::is_enum_v<R> || std::is_enum_v<C>))304
constexpr const T& operator()(R row, C col) const {305
return (*this)(static_cast<std::size_t>(::cheatah::ndarray::subscript_index(row)),306
static_cast<std::size_t>(::cheatah::ndarray::subscript_index(col)));307
}309
/**310
* A pointer to the elements, contiguous — column-major for a matrix, which is exactly the order a311
* GPU uniform, a push constant or a BLAS call expects, so an upload is a copy not a transpose.312
* @return the first element's address.313
* @complexity O(1).314
* @alloc none.315
* @test Fixarray.Data316
*/317
constexpr T* data() { return data_.data(); }319
/**320
* A pointer to the elements, contiguous and column-major for a matrix (read-only).321
* @return the first element's address.322
* @complexity O(1).323
* @alloc none.324
* @test Fixarray.Data325
*/326
constexpr const T* data() const { return data_.data(); }328
/**329
* Elementwise equality. Exact, as `==` on the elements is exact — no tolerance is applied to330
* floating-point values.331
* @param other the array to compare with.332
* @return true iff every element matches.333
* @complexity O(size).334
* @alloc none.335
* @test Fixarray.Equality336
*/337
constexpr bool operator==(const Fixed& other) const = default;339
/**340
* Add @p other elementwise, in place.341
* @param other the array to add.342
* @return a reference to this array.343
* @complexity O(size).344
* @alloc none.345
* @test Fixarray.Arithmetic346
*/347
constexpr Fixed& operator+=(const Fixed& other) {348
for (std::size_t i = 0; i < size; ++i) { data_[i] += other.data_[i]; }349
return *this;350
}352
/**353
* Subtract @p other elementwise, in place.354
* @param other the array to subtract.355
* @return a reference to this array.356
* @complexity O(size).357
* @alloc none.358
* @test Fixarray.Arithmetic359
*/360
constexpr Fixed& operator-=(const Fixed& other) {361
for (std::size_t i = 0; i < size; ++i) { data_[i] -= other.data_[i]; }362
return *this;363
}365
/**366
* Scale every element by @p scalar, in place.367
* @param scalar the factor.368
* @return a reference to this array.369
* @complexity O(size).370
* @alloc none.371
* @test Fixarray.Arithmetic372
*/373
constexpr Fixed& operator*=(T scalar) {374
for (std::size_t i = 0; i < size; ++i) { data_[i] *= scalar; }375
return *this;376
}378
/**379
* Divide every element by @p scalar, in place.380
* @param scalar the divisor.381
* @return a reference to this array.382
* @complexity O(size).383
* @alloc none.384
* @test Fixarray.Arithmetic385
*/386
constexpr Fixed& operator/=(T scalar) {387
for (std::size_t i = 0; i < size; ++i) { data_[i] /= scalar; }388
return *this;389
}391
/**392
* Elementwise sum.393
* @param a,b the arrays to add.394
* @return `a + b`.395
* @complexity O(size).396
* @alloc none.397
* @test Fixarray.Arithmetic398
*/399
friend constexpr Fixed operator+(Fixed a, const Fixed& b) { return a += b; }401
/**402
* Elementwise difference.403
* @param a,b the arrays to subtract.404
* @return `a - b`.405
* @complexity O(size).406
* @alloc none.407
* @test Fixarray.Arithmetic408
*/409
friend constexpr Fixed operator-(Fixed a, const Fixed& b) { return a -= b; }411
/**412
* Negation.413
* @param a the array to negate.414
* @return `-a`.415
* @complexity O(size).416
* @alloc none.417
* @test Fixarray.Arithmetic418
*/419
friend constexpr Fixed operator-(Fixed a) { return a *= static_cast<T>(-1); }421
/**422
* Scale by a scalar.423
* @param a the array. @param scalar the factor.424
* @return `a * scalar`.425
* @complexity O(size).426
* @alloc none.427
* @test Fixarray.Arithmetic428
*/429
friend constexpr Fixed operator*(Fixed a, T scalar) { return a *= scalar; }431
/**432
* Scale by a scalar.433
* @param scalar the factor. @param a the array.434
* @return `scalar * a`.435
* @complexity O(size).436
* @alloc none.437
* @test Fixarray.Arithmetic438
*/439
friend constexpr Fixed operator*(T scalar, Fixed a) { return a *= scalar; }441
/**442
* Divide by a scalar.443
* @param a the array. @param scalar the divisor.444
* @return `a / scalar`.445
* @complexity O(size).446
* @alloc none.447
* @test Fixarray.Arithmetic448
*/449
friend constexpr Fixed operator/(Fixed a, T scalar) { return a /= scalar; }451
private:452
/// Adopt an already-built element buffer, skipping the zero-initialization of the default453
/// constructor. Private: the buffer's layout (column-major for a matrix) is an implementation454
/// detail that only the members below may rely on.455
explicit constexpr Fixed(const std::array<T, size>& elements) : data_(elements) {}457
/// @ref identity's worker: emits each element exactly once, with no zeroing pass.458
/// @tparam I the flat indices 0 … size-1.459
/// @return the identity matrix.460
template <std::size_t... I>461
static constexpr Fixed identity_impl(std::index_sequence<I...> /*unused*/)462
requires(rank == 2 && rows == cols)463
{464
return Fixed(std::array<T, size>{(I % (rows + 1) == 0 ? T{1} : T{0})...});465
}467
/// @ref from_indices's worker: aggregate-initialises the buffer from `f(0) … f(size-1)`, fully468
/// unrolled, so there is no loop, no zeroing, and no pointer through which the operands alias.469
/// @tparam F the element-producing callable.470
/// @tparam I the flat indices 0 … size-1.471
/// @param f produces each element from its flat index.472
/// @return the array of `f(i)`.473
template <class F, std::size_t... I>474
static constexpr Fixed from_indices_impl(F&& f, std::index_sequence<I...> /*unused*/) { // NOLINT(cppcoreguidelines-missing-std-forward): f is invoked once per index; forwarding inside the pack expansion would move it repeatedly475
return Fixed(std::array<T, size>{static_cast<T>(f(I))...});476
}478
/// The elements, inline: a vector in order, a matrix column by column. Zero by default.479
std::array<T, size> data_{};480
};482
/// A fixed-extent vector of @p N elements.483
/// @tparam T the element type. @tparam N the length.484
template <ndarray::Field T, std::size_t N>485
using Vec = Fixed<T, N>;487
/// A fixed-extent matrix of @p R rows and @p C columns — column-major storage, mathematical488
/// `(row, col)` indexing (see the file doc).489
/// @tparam T the element type. @tparam R the rows. @tparam C the columns.490
template <ndarray::Field T, std::size_t R, std::size_t C>491
using Mat = Fixed<T, R, C>;493
/// A 2-D vector of `float`.494
using vec2f = Vec<float, 2>;495
/// A 3-D vector of `float` — a direction, a position, a colour.496
using vec3f = Vec<float, 3>;497
/// A 4-D vector of `float` — a homogeneous point, an RGBA colour.498
using vec4f = Vec<float, 4>;499
/// A 2-D vector of `double`.500
using vec2d = Vec<double, 2>;501
/// A 3-D vector of `double`.502
using vec3d = Vec<double, 3>;503
/// A 4-D vector of `double`.504
using vec4d = Vec<double, 4>;506
/// A 2×2 matrix of `float`.507
using mat2f = Mat<float, 2, 2>;508
/// A 3×3 matrix of `float` — a rotation, or a normal matrix.509
using mat3f = Mat<float, 3, 3>;510
/// A 4×4 matrix of `float` — a transform; exactly the 64 bytes of a push constant.511
using mat4f = Mat<float, 4, 4>;512
/// A 2×2 matrix of `double`.513
using mat2d = Mat<double, 2, 2>;514
/// A 3×3 matrix of `double`.515
using mat3d = Mat<double, 3, 3>;516
/// A 4×4 matrix of `double`.517
using mat4d = Mat<double, 4, 4>;519
namespace detail {521
/**522
* Sum @p n elements PAIRWISE rather than left to right. A serial `sum += x[i]` chains each add on523
* the previous one, so the loop runs at the latency of an addition; halving the array instead lets524
* independent adds issue together, and — the reason numerics people reach for it — the rounding525
* error grows as O(log n) instead of O(n).526
* @tparam T the element type. @tparam N the array length.527
* @param values the products to sum.528
* @return their sum.529
* @complexity O(N).530
* @alloc none.531
* @test Fixarray.DotAndCross532
*/533
template <ndarray::Field T, std::size_t N>534
constexpr T pairwise_sum(const std::array<T, N>& values) {535
if constexpr (N == 1) {536
return values[0];537
} else if constexpr (N == 2) {538
return values[0] + values[1];539
} else if constexpr (N == 3) {540
return (values[0] + values[1]) + values[2];541
} else if constexpr (N == 4) {542
return (values[0] + values[1]) + (values[2] + values[3]);543
} else {544
constexpr std::size_t half = N / 2;545
std::array<T, half> lo{};546
std::array<T, N - half> hi{};547
for (std::size_t i = 0; i < half; ++i) { lo[i] = values[i]; }548
for (std::size_t i = half; i < N; ++i) { hi[i - half] = values[i]; }549
return pairwise_sum(lo) + pairwise_sum(hi);550
}551
}553
} // namespace detail555
/**556
* Inner product of two vectors — Σ aᵢbᵢ, the same quantity @ref dot(const NDArray&, const NDArray&)557
* computes, without the allocation. Summed pairwise, so it is both faster and more accurate than a558
* left-to-right accumulation.559
* @tparam T the element type. @tparam N the length.560
* @param a,b the vectors.561
* @return the inner product.562
* @complexity O(N).563
* @alloc none.564
* @test Fixarray.DotAndCross565
*/566
template <ndarray::Field T, std::size_t N>567
constexpr T dot(const Vec<T, N>& a, const Vec<T, N>& b) {568
// Two regimes, split by what the compiler does with the products. Below 4, an odd width — a569
// 3-vector of doubles is 24 bytes — spills a products array to the stack, so the terms are570
// written out and stay in registers. At 4 and above the array is a whole SIMD register (or a571
// clean multiple), and the loop packs the products into one `mulps`/`mulpd` instead of N scalar572
// multiplies — which is how this *beats* a scalar dot rather than merely matching it. Both paths573
// sum pairwise, so the associativity (and the rounding) is identical either way.574
if constexpr (N == 1) {575
return a[0] * b[0];576
} else if constexpr (N == 2) {577
return a[0] * b[0] + a[1] * b[1];578
} else if constexpr (N == 3) {579
return (a[0] * b[0] + a[1] * b[1]) + a[2] * b[2];580
} else {581
std::array<T, N> products{};582
for (std::size_t i = 0; i < N; ++i) { products[i] = a[i] * b[i]; }583
return detail::pairwise_sum(products);584
}585
}587
/**588
* Cross product of two 3-vectors — the vector perpendicular to both, right-handed.589
* @tparam T the element type.590
* @param a,b the vectors.591
* @return `a × b`.592
* @complexity O(1).593
* @alloc none.594
* @test Fixarray.DotAndCross595
*/596
template <ndarray::Field T>597
constexpr Vec<T, 3> cross(const Vec<T, 3>& a, const Vec<T, 3>& b) {598
return Vec<T, 3>{a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],599
a[0] * b[1] - a[1] * b[0]};600
}602
/**603
* The squared Euclidean length of a vector — `dot(v, v)`. Prefer it to @ref norm when only comparing604
* lengths: it avoids the square root.605
* @tparam T the element type. @tparam N the length.606
* @param v the vector.607
* @return `Σ vᵢ²`.608
* @complexity O(N).609
* @alloc none.610
* @warning for a complex element this is the bilinear `Σ vᵢ²` (matching @ref dot), not611
* the Hermitian `Σ |vᵢ|²` — it is the squared Euclidean LENGTH only for real612
* elements.613
* @test Fixarray.NormAndNormalize614
*/615
template <ndarray::Field T, std::size_t N>616
constexpr T squared_norm(const Vec<T, N>& v) {617
return dot(v, v);618
}620
/**621
* The Euclidean length of a vector.622
* @tparam T the element type; floating-point, since the result is a root.623
* @tparam N the length.624
* @param v the vector.625
* @return `sqrt(Σ vᵢ²)`.626
* @complexity O(N).627
* @alloc none.628
* @test Fixarray.NormAndNormalize629
*/630
template <ndarray::FloatingPoint T, std::size_t N>631
T norm(const Vec<T, N>& v) {632
return std::sqrt(squared_norm(v));633
}635
/**636
* The unit vector pointing the same way as @p v.637
* @tparam T the element type; floating-point.638
* @tparam N the length.639
* @param v the vector; must not be the zero vector.640
* @return `v / norm(v)`.641
* @throws std::domain_error when @p v has zero length, since it has no direction.642
* @complexity O(N).643
* @alloc none.644
* @test Fixarray.NormAndNormalize645
*/646
template <ndarray::FloatingPoint T, std::size_t N>647
Vec<T, N> normalize(const Vec<T, N>& v) {648
const T squared = squared_norm(v);649
if (squared == T{0}) { throw std::domain_error("fixarray::normalize: the zero vector has no direction"); }650
// One reciprocal, then N multiplies. Dividing each component instead costs N divides, and a651
// divide is roughly three times the latency of a multiply.652
const T inverse_length = T{1} / std::sqrt(squared);653
return v * inverse_length;654
}656
/**657
* The transpose of a matrix — rows become columns.658
* @tparam T the element type. @tparam R the rows. @tparam C the columns.659
* @param m the matrix.660
* @return the `C×R` transpose.661
* @complexity O(R·C).662
* @alloc none.663
* @test Fixarray.TransposeAndTrace664
*/665
template <ndarray::Field T, std::size_t R, std::size_t C>666
constexpr Mat<T, C, R> transpose(const Mat<T, R, C>& m) {667
Mat<T, C, R> result;668
for (std::size_t r = 0; r < R; ++r) {669
for (std::size_t c = 0; c < C; ++c) { result(c, r) = m(r, c); }670
}671
return result;672
}674
/**675
* The trace of a square matrix — the sum of its diagonal.676
* @tparam T the element type. @tparam N the dimension.677
* @param m the matrix.678
* @return `Σ mᵢᵢ`.679
* @complexity O(N).680
* @alloc none.681
* @test Fixarray.TransposeAndTrace682
*/683
template <ndarray::Field T, std::size_t N>684
constexpr T trace(const Mat<T, N, N>& m) {685
T sum{};686
for (std::size_t i = 0; i < N; ++i) { sum += m(i, i); }687
return sum;688
}690
/**691
* Matrix product — the same `A·B` @ref matmul(const NDArray&, const NDArray&) computes, with the692
* shapes checked by the compiler rather than at runtime.693
* @tparam T the element type. @tparam R the rows of @p a. @tparam K the shared dimension.694
* @tparam C the columns of @p b.695
* @param a,b the matrices.696
* @return the `R×C` product.697
* @complexity O(R·K·C).698
* @alloc none.699
* @test Fixarray.Matmul700
*/701
template <ndarray::Field T, std::size_t R, std::size_t K, std::size_t C>702
constexpr Mat<T, R, C> matmul(const Mat<T, R, K>& a, const Mat<T, K, C>& b) {703
Mat<T, R, C> result;704
// Column c of the product is Σₖ (column k of a) · b(k, c): a sum of SCALED COLUMNS. With705
// column-major storage each of those columns is contiguous, so the inner loop is a plain vertical706
// multiply-add that vectorizes. The first term SEEDS the column rather than adding to a zeroed707
// one, which spares a store-then-reload of the accumulator.708
for (std::size_t c = 0; c < C; ++c) {709
const T first = b(0, c);710
for (std::size_t r = 0; r < R; ++r) { result(r, c) = first * a(r, 0); }711
for (std::size_t k = 1; k < K; ++k) {712
const T scale = b(k, c);713
for (std::size_t r = 0; r < R; ++r) { result(r, c) += scale * a(r, k); }714
}715
}716
return result;717
}719
/**720
* Matrix product, spelled `a * b`.721
* @tparam T the element type. @tparam R the rows of @p a. @tparam K the shared dimension.722
* @tparam C the columns of @p b.723
* @param a,b the matrices.724
* @return `matmul(a, b)`.725
* @complexity O(R·K·C).726
* @alloc none.727
* @test Fixarray.Matmul728
*/729
template <ndarray::Field T, std::size_t R, std::size_t K, std::size_t C>730
constexpr Mat<T, R, C> operator*(const Mat<T, R, K>& a, const Mat<T, K, C>& b) {731
return matmul(a, b);732
}734
/**735
* Transform a vector by a matrix — `A·v`, treating @p v as a column.736
* @tparam T the element type. @tparam R the rows. @tparam C the columns, and @p v's length.737
* @param m the matrix. @param v the vector.738
* @return the `R`-vector `m · v`.739
* @complexity O(R·C).740
* @alloc none.741
* @test Fixarray.Matmul742
*/743
template <ndarray::Field T, std::size_t R, std::size_t C>744
constexpr Vec<T, R> operator*(const Mat<T, R, C>& m, const Vec<T, C>& v) {745
// `m · v` is Σⱼ (column j of m) · v[j] — a sum of scaled columns, not a stack of row dot746
// products. Column-major storage makes each column contiguous, so this is a vertical747
// multiply-add with no horizontal reduction and no shuffles. The first column seeds the result748
// rather than adding into a zeroed one, which spares a pass over it.749
Vec<T, R> result;750
const T first = v[0];751
for (std::size_t r = 0; r < R; ++r) { result[r] = first * m(r, 0); }752
for (std::size_t c = 1; c < C; ++c) {753
const T scale = v[c];754
for (std::size_t r = 0; r < R; ++r) { result[r] += scale * m(r, c); }755
}756
return result;757
}759
/**760
* The determinant of a 2×2 matrix, in closed form.761
* @tparam T the element type.762
* @param m the matrix.763
* @return `det(m)`.764
* @complexity O(1).765
* @alloc none.766
* @test Fixarray.DeterminantAndInverse767
*/768
template <ndarray::Field T>769
constexpr T determinant(const Mat<T, 2, 2>& m) {770
return m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0);771
}773
/**774
* The determinant of a 3×3 matrix, by the rule of Sarrus.775
* @tparam T the element type.776
* @param m the matrix.777
* @return `det(m)`.778
* @complexity O(1).779
* @alloc none.780
* @test Fixarray.DeterminantAndInverse781
*/782
template <ndarray::Field T>783
constexpr T determinant(const Mat<T, 3, 3>& m) {784
return m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) -785
m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0)) +786
m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));787
}789
/**790
* The determinant of a 4×4 matrix, by cofactor expansion on 2×2 minors — the form a transform791
* matrix meets, and cheaper than an LU factorization at this size.792
* @tparam T the element type.793
* @param m the matrix.794
* @return `det(m)`.795
* @complexity O(1).796
* @alloc none.797
* @test Fixarray.DeterminantAndInverse798
*/799
template <ndarray::Field T>800
constexpr T determinant(const Mat<T, 4, 4>& m) {801
const T s0 = m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1);802
const T s1 = m(0, 0) * m(1, 2) - m(1, 0) * m(0, 2);803
const T s2 = m(0, 0) * m(1, 3) - m(1, 0) * m(0, 3);804
const T s3 = m(0, 1) * m(1, 2) - m(1, 1) * m(0, 2);805
const T s4 = m(0, 1) * m(1, 3) - m(1, 1) * m(0, 3);806
const T s5 = m(0, 2) * m(1, 3) - m(1, 2) * m(0, 3);808
const T c5 = m(2, 2) * m(3, 3) - m(3, 2) * m(2, 3);809
const T c4 = m(2, 1) * m(3, 3) - m(3, 1) * m(2, 3);810
const T c3 = m(2, 1) * m(3, 2) - m(3, 1) * m(2, 2);811
const T c2 = m(2, 0) * m(3, 3) - m(3, 0) * m(2, 3);812
const T c1 = m(2, 0) * m(3, 2) - m(3, 0) * m(2, 2);813
const T c0 = m(2, 0) * m(3, 1) - m(3, 0) * m(2, 1);815
return s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0;816
}818
/**819
* The inverse of a 2×2 matrix, in closed form.820
* @tparam T the element type; floating-point, since the inverse divides.821
* @param m the matrix; must be non-singular.822
* @return `m⁻¹`.823
* @throws std::domain_error when @p m is singular (zero determinant).824
* @complexity O(1).825
* @alloc none.826
* @test Fixarray.DeterminantAndInverse827
*/828
template <ndarray::FloatingPoint T>829
constexpr Mat<T, 2, 2> inverse(const Mat<T, 2, 2>& m) {830
const T det = determinant(m);831
if (det == T{0}) { throw std::domain_error("fixarray::inverse: the matrix is singular"); }832
const T inv_det = T{1} / det;833
Mat<T, 2, 2> result;834
result(0, 0) = m(1, 1) * inv_det;835
result(0, 1) = -m(0, 1) * inv_det;836
result(1, 0) = -m(1, 0) * inv_det;837
result(1, 1) = m(0, 0) * inv_det;838
return result;839
}841
/**842
* The inverse of a 3×3 matrix, by its adjugate — the normal matrix a renderer needs.843
* @tparam T the element type; floating-point.844
* @param m the matrix; must be non-singular.845
* @return `m⁻¹`.846
* @throws std::domain_error when @p m is singular (zero determinant).847
* @complexity O(1).848
* @alloc none.849
* @test Fixarray.DeterminantAndInverse850
*/851
template <ndarray::FloatingPoint T>852
constexpr Mat<T, 3, 3> inverse(const Mat<T, 3, 3>& m) {853
const T det = determinant(m);854
if (det == T{0}) { throw std::domain_error("fixarray::inverse: the matrix is singular"); }855
const T inv_det = T{1} / det;856
Mat<T, 3, 3> result;857
result(0, 0) = (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) * inv_det;858
result(0, 1) = (m(0, 2) * m(2, 1) - m(0, 1) * m(2, 2)) * inv_det;859
result(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) * inv_det;860
result(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) * inv_det;861
result(1, 1) = (m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0)) * inv_det;862
result(1, 2) = (m(0, 2) * m(1, 0) - m(0, 0) * m(1, 2)) * inv_det;863
result(2, 0) = (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)) * inv_det;864
result(2, 1) = (m(0, 1) * m(2, 0) - m(0, 0) * m(2, 1)) * inv_det;865
result(2, 2) = (m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0)) * inv_det;866
return result;867
}869
/**870
* The inverse of a 4×4 matrix, by its adjugate over the 2×2 minors — the transform a camera871
* inverts every frame.872
* @tparam T the element type; floating-point.873
* @param m the matrix; must be non-singular.874
* @return `m⁻¹`.875
* @throws std::domain_error when @p m is singular (zero determinant).876
* @complexity O(1).877
* @alloc none.878
* @test Fixarray.DeterminantAndInverse879
*/880
template <ndarray::FloatingPoint T>881
constexpr Mat<T, 4, 4> inverse(const Mat<T, 4, 4>& m) {882
const T s0 = m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1);883
const T s1 = m(0, 0) * m(1, 2) - m(1, 0) * m(0, 2);884
const T s2 = m(0, 0) * m(1, 3) - m(1, 0) * m(0, 3);885
const T s3 = m(0, 1) * m(1, 2) - m(1, 1) * m(0, 2);886
const T s4 = m(0, 1) * m(1, 3) - m(1, 1) * m(0, 3);887
const T s5 = m(0, 2) * m(1, 3) - m(1, 2) * m(0, 3);889
const T c5 = m(2, 2) * m(3, 3) - m(3, 2) * m(2, 3);890
const T c4 = m(2, 1) * m(3, 3) - m(3, 1) * m(2, 3);891
const T c3 = m(2, 1) * m(3, 2) - m(3, 1) * m(2, 2);892
const T c2 = m(2, 0) * m(3, 3) - m(3, 0) * m(2, 3);893
const T c1 = m(2, 0) * m(3, 2) - m(3, 0) * m(2, 2);894
const T c0 = m(2, 0) * m(3, 1) - m(3, 0) * m(2, 1);896
const T det = s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0;897
if (det == T{0}) { throw std::domain_error("fixarray::inverse: the matrix is singular"); }898
const T d = T{1} / det;900
Mat<T, 4, 4> r;901
r(0, 0) = (m(1, 1) * c5 - m(1, 2) * c4 + m(1, 3) * c3) * d;902
r(0, 1) = (-m(0, 1) * c5 + m(0, 2) * c4 - m(0, 3) * c3) * d;903
r(0, 2) = (m(3, 1) * s5 - m(3, 2) * s4 + m(3, 3) * s3) * d;904
r(0, 3) = (-m(2, 1) * s5 + m(2, 2) * s4 - m(2, 3) * s3) * d;906
r(1, 0) = (-m(1, 0) * c5 + m(1, 2) * c2 - m(1, 3) * c1) * d;907
r(1, 1) = (m(0, 0) * c5 - m(0, 2) * c2 + m(0, 3) * c1) * d;908
r(1, 2) = (-m(3, 0) * s5 + m(3, 2) * s2 - m(3, 3) * s1) * d;909
r(1, 3) = (m(2, 0) * s5 - m(2, 2) * s2 + m(2, 3) * s1) * d;911
r(2, 0) = (m(1, 0) * c4 - m(1, 1) * c2 + m(1, 3) * c0) * d;912
r(2, 1) = (-m(0, 0) * c4 + m(0, 1) * c2 - m(0, 3) * c0) * d;913
r(2, 2) = (m(3, 0) * s4 - m(3, 1) * s2 + m(3, 3) * s0) * d;914
r(2, 3) = (-m(2, 0) * s4 + m(2, 1) * s2 - m(2, 3) * s0) * d;916
r(3, 0) = (-m(1, 0) * c3 + m(1, 1) * c1 - m(1, 2) * c0) * d;917
r(3, 1) = (m(0, 0) * c3 - m(0, 1) * c1 + m(0, 2) * c0) * d;918
r(3, 2) = (-m(3, 0) * s3 + m(3, 1) * s1 - m(3, 2) * s0) * d;919
r(3, 3) = (m(2, 0) * s3 - m(2, 1) * s1 + m(2, 2) * s0) * d;920
return r;921
}923
// ---- Geometry: the operations a renderer and a physics solver reach for ------------------------924
// These are the GLSL/GLM geometric builtins, by their standard names, over @ref Fixed vectors: the925
// same mathematics, evaluated in registers with no allocation. They reuse the products above, so a926
// change to @ref dot or the operators reaches them too.928
/**929
* The Euclidean distance between two points — `norm(a - b)`.930
* @tparam T the element type; floating-point, since the result is a root.931
* @tparam N the dimension.932
* @param a,b the points.933
* @return `‖a − b‖`.934
* @complexity O(N). @alloc none.935
* @test Fixarray.Geometry936
*/937
template <ndarray::FloatingPoint T, std::size_t N>938
T distance(const Vec<T, N>& a, const Vec<T, N>& b) {939
return norm(a - b);940
}942
/**943
* The squared distance between two points — `squared_norm(a - b)`. Prefer it to @ref distance when944
* only comparing distances: it skips the square root.945
* @tparam T the element type. @tparam N the dimension.946
* @param a,b the points.947
* @return `‖a − b‖²`.948
* @complexity O(N). @alloc none.949
* @test Fixarray.Geometry950
*/951
template <ndarray::Field T, std::size_t N>952
constexpr T distance_squared(const Vec<T, N>& a, const Vec<T, N>& b) {953
return squared_norm(a - b);954
}956
/**957
* Reflect an incident vector about a surface normal — `I − 2 (N·I) N`, the GLSL `reflect`. @p normal958
* is assumed unit length, as GLSL requires.959
* @tparam T the element type. @tparam N the dimension.960
* @param incident the incoming vector.961
* @param normal the unit surface normal.962
* @return the reflected vector.963
* @complexity O(N). @alloc none.964
* @test Fixarray.Geometry965
*/966
template <ndarray::Field T, std::size_t N>967
constexpr Vec<T, N> reflect(const Vec<T, N>& incident, const Vec<T, N>& normal) {968
return incident - (T{2} * dot(normal, incident)) * normal;969
}971
/**972
* Refract an incident vector through a surface — the GLSL `refract`. @p incident and @p normal are973
* assumed unit length. On total internal reflection (a negative radicand) the result is the zero974
* vector, exactly as GLSL specifies.975
* @tparam T the element type; floating-point.976
* @tparam N the dimension.977
* @param incident the unit incoming vector.978
* @param normal the unit surface normal.979
* @param eta the ratio of indices of refraction (source over destination).980
* @return the refracted vector, or the zero vector under total internal reflection.981
* @complexity O(N). @alloc none.982
* @test Fixarray.Geometry983
*/984
template <ndarray::FloatingPoint T, std::size_t N>985
Vec<T, N> refract(const Vec<T, N>& incident, const Vec<T, N>& normal, T eta) {986
const T cos_i = dot(normal, incident);987
const T k = T{1} - eta * eta * (T{1} - cos_i * cos_i);988
if (k < T{0}) { return Vec<T, N>{}; }989
return eta * incident - (eta * cos_i + std::sqrt(k)) * normal;990
}992
/**993
* Orient a normal to face a viewer — the GLSL `faceforward`: return @p n when @p reference points994
* against the incident direction (`dot(reference, incident) < 0`), `-n` otherwise. Used to keep a995
* surface normal on the camera's side.996
* @tparam T the element type. @tparam N the dimension.997
* @param n the normal to orient.998
* @param incident the incident vector.999
* @param reference the reference normal the result is oriented against.1000
* @return @p n or `-n`.1001
* @complexity O(N). @alloc none.1002
* @test Fixarray.Geometry1003
*/1004
template <ndarray::Numeric T, std::size_t N>1005
constexpr Vec<T, N> faceforward(const Vec<T, N>& n, const Vec<T, N>& incident,1006
const Vec<T, N>& reference) {1007
return dot(reference, incident) < T{0} ? n : -n;1008
}1010
// ---- Component-wise functions: the GLSL/GLM "common" builtins over a whole array ----------------1011
// Each applies elementwise to every element of a @ref Fixed — a vector or a matrix alike — so a1012
// renderer clamps a colour, a physics step limits a velocity, and a noise field mixes two samples in1013
// the same vocabulary. They read and write the flat buffer, so they are correct whatever the storage1014
// order, and they copy their first argument rather than zero a result and overwrite it.1016
/**1017
* The absolute value of every element.1018
* @tparam T the element type; a real number, so `< 0` is meaningful.1019
* @tparam Dims the extents.1020
* @param x the array.1021
* @return `|x|` elementwise.1022
* @complexity O(size). @alloc none.1023
* @test Fixarray.CommonUnary1024
*/1025
template <ndarray::Numeric T, std::size_t... Dims>1026
constexpr Fixed<T, Dims...> abs(const Fixed<T, Dims...>& x) {1027
return Fixed<T, Dims...>::from_indices([&](std::size_t i) {1028
const T v = x.data()[i];1029
return v < T{0} ? -v : v;1030
});1031
}1033
/**1034
* The sign of every element: `-1`, `0`, or `+1`.1035
* @tparam T the element type; a real number.1036
* @tparam Dims the extents.1037
* @param x the array.1038
* @return the elementwise sign.1039
* @complexity O(size). @alloc none.1040
* @test Fixarray.CommonUnary1041
*/1042
template <ndarray::Numeric T, std::size_t... Dims>1043
constexpr Fixed<T, Dims...> sign(const Fixed<T, Dims...>& x) {1044
return Fixed<T, Dims...>::from_indices([&](std::size_t i) {1045
const T v = x.data()[i];1046
return static_cast<T>((T{0} < v) - (v < T{0}));1047
});1048
}1050
/**1051
* The smaller of each corresponding pair of elements.1052
* @tparam T the element type; a real number.1053
* @tparam Dims the extents.1054
* @param a,b the arrays.1055
* @return `min(aᵢ, bᵢ)` elementwise.1056
* @complexity O(size). @alloc none.1057
* @test Fixarray.MinMaxClamp1058
*/1059
template <ndarray::Numeric T, std::size_t... Dims>1060
constexpr Fixed<T, Dims...> min(const Fixed<T, Dims...>& a, const Fixed<T, Dims...>& b) {1061
return Fixed<T, Dims...>::from_indices([&](std::size_t i) {1062
const T ai = a.data()[i];1063
const T bi = b.data()[i];1064
return ai < bi ? ai : bi; // a branchless min lowers to minps/minpd1065
});1066
}1068
/**1069
* Each element capped at the scalar @p s — `min(xᵢ, s)`.1070
* @tparam T the element type; a real number.1071
* @tparam Dims the extents.1072
* @param x the array. @param s the ceiling applied to every element.1073
* @return `min(xᵢ, s)` elementwise.1074
* @complexity O(size). @alloc none.1075
* @test Fixarray.MinMaxClamp1076
*/1077
template <ndarray::Numeric T, std::size_t... Dims>1078
constexpr Fixed<T, Dims...> min(const Fixed<T, Dims...>& x, T s) {1079
return Fixed<T, Dims...>::from_indices([&](std::size_t i) {1080
const T v = x.data()[i];1081
return v < s ? v : s;1082
});1083
}1085
/**1086
* The larger of each corresponding pair of elements.1087
* @tparam T the element type; a real number.1088
* @tparam Dims the extents.1089
* @param a,b the arrays.1090
* @return `max(aᵢ, bᵢ)` elementwise.1091
* @complexity O(size). @alloc none.1092
* @test Fixarray.MinMaxClamp1093
*/1094
template <ndarray::Numeric T, std::size_t... Dims>1095
constexpr Fixed<T, Dims...> max(const Fixed<T, Dims...>& a, const Fixed<T, Dims...>& b) {1096
return Fixed<T, Dims...>::from_indices([&](std::size_t i) {1097
const T ai = a.data()[i];1098
const T bi = b.data()[i];1099
return ai < bi ? bi : ai; // branchless max -> maxps/maxpd1100
});1101
}1103
/**1104
* Each element raised to the scalar @p s — `max(xᵢ, s)`.1105
* @tparam T the element type; a real number.1106
* @tparam Dims the extents.1107
* @param x the array. @param s the floor applied to every element.1108
* @return `max(xᵢ, s)` elementwise.1109
* @complexity O(size). @alloc none.1110
* @test Fixarray.MinMaxClamp1111
*/1112
template <ndarray::Numeric T, std::size_t... Dims>1113
constexpr Fixed<T, Dims...> max(const Fixed<T, Dims...>& x, T s) {1114
return Fixed<T, Dims...>::from_indices([&](std::size_t i) {1115
const T v = x.data()[i];1116
return v < s ? s : v;1117
});1118
}1120
/**1121
* Constrain every element to `[lo, hi]` — the GLSL `clamp` with scalar bounds, the common case of1122
* pinning a colour to `[0, 1]`.1123
* @tparam T the element type; a real number.1124
* @tparam Dims the extents.1125
* @param x the array. @param lo the lower bound. @param hi the upper bound.1126
* @return `min(max(xᵢ, lo), hi)` elementwise.1127
* @complexity O(size). @alloc none.1128
* @test Fixarray.MinMaxClamp1129
*/1130
template <ndarray::Numeric T, std::size_t... Dims>1131
constexpr Fixed<T, Dims...> clamp(const Fixed<T, Dims...>& x, T lo, T hi) {1132
return Fixed<T, Dims...>::from_indices([&](std::size_t i) {1133
const T v = x.data()[i];1134
const T low = v < lo ? lo : v;1135
return hi < low ? hi : low; // min(max(v, lo), hi), branchless1136
});1137
}1139
/**1140
* Constrain every element between the corresponding bounds — the GLSL `clamp` with per-element1141
* bounds.1142
* @tparam T the element type; a real number.1143
* @tparam Dims the extents.1144
* @param x the array. @param lo the lower bounds. @param hi the upper bounds.1145
* @return `min(max(xᵢ, loᵢ), hiᵢ)` elementwise.1146
* @complexity O(size). @alloc none.1147
* @test Fixarray.MinMaxClamp1148
*/1149
template <ndarray::Numeric T, std::size_t... Dims>1150
constexpr Fixed<T, Dims...> clamp(const Fixed<T, Dims...>& x, const Fixed<T, Dims...>& lo,1151
const Fixed<T, Dims...>& hi) {1152
return Fixed<T, Dims...>::from_indices([&](std::size_t i) {1153
const T v = x.data()[i];1154
const T l = lo.data()[i];1155
const T h = hi.data()[i];1156
const T low = v < l ? l : v;1157
return h < low ? h : low;1158
});1159
}1161
/**1162
* Linear interpolation — the GLSL `mix`: `a (1 − t) + b t`, with a scalar blend @p t (0 gives @p a,1163
* 1 gives @p b). Composed from the operators, so it inherits their vectorization.1164
* @tparam T the element type; floating-point.1165
* @tparam Dims the extents.1166
* @param a,b the endpoints. @param t the blend factor.1167
* @return the interpolated array.1168
* @complexity O(size). @alloc none.1169
* @test Fixarray.MixStep1170
*/1171
template <ndarray::FloatingPoint T, std::size_t... Dims>1172
constexpr Fixed<T, Dims...> mix(const Fixed<T, Dims...>& a, const Fixed<T, Dims...>& b, T t) {1173
return a * (T{1} - t) + b * t;1174
}1176
/**1177
* Linear interpolation with a per-element blend — the GLSL `mix` whose factor @p t is an array.1178
* @tparam T the element type; floating-point.1179
* @tparam Dims the extents.1180
* @param a,b the endpoints. @param t the per-element blend factors.1181
* @return the interpolated array.1182
* @complexity O(size). @alloc none.1183
* @test Fixarray.MixStep1184
*/1185
template <ndarray::FloatingPoint T, std::size_t... Dims>1186
constexpr Fixed<T, Dims...> mix(const Fixed<T, Dims...>& a, const Fixed<T, Dims...>& b,1187
const Fixed<T, Dims...>& t) {1188
return Fixed<T, Dims...>::from_indices(1189
[&](std::size_t i) { return a.data()[i] * (T{1} - t.data()[i]) + b.data()[i] * t.data()[i]; });1190
}1192
/**1193
* A step at @p edge — the GLSL `step`: `0` where an element is below @p edge, `1` at or above.1194
* @tparam T the element type; a real number.1195
* @tparam Dims the extents.1196
* @param edge the threshold. @param x the array.1197
* @return `xᵢ < edge ? 0 : 1` elementwise.1198
* @complexity O(size). @alloc none.1199
* @test Fixarray.MixStep1200
*/1201
template <ndarray::Numeric T, std::size_t... Dims>1202
constexpr Fixed<T, Dims...> step(T edge, const Fixed<T, Dims...>& x) {1203
return Fixed<T, Dims...>::from_indices(1204
[&](std::size_t i) { return x.data()[i] < edge ? T{0} : T{1}; });1205
}1207
/**1208
* A smooth Hermite transition from 0 to 1 across `[edge0, edge1]` — the GLSL `smoothstep`, with1209
* everything below @p edge0 giving 0 and everything above @p edge1 giving 1.1210
* @tparam T the element type; floating-point.1211
* @tparam Dims the extents.1212
* @param edge0 the lower edge. @param edge1 the upper edge. @param x the array.1213
* @return the smoothstepped array.1214
* @complexity O(size). @alloc none.1215
* @test Fixarray.MixStep1216
*/1217
template <ndarray::FloatingPoint T, std::size_t... Dims>1218
constexpr Fixed<T, Dims...> smoothstep(T edge0, T edge1, const Fixed<T, Dims...>& x) {1219
return Fixed<T, Dims...>::from_indices([&](std::size_t i) {1220
T t = (x.data()[i] - edge0) / (edge1 - edge0);1221
t = t < T{0} ? T{0} : (T{1} < t ? T{1} : t); // NOLINT(readability-avoid-nested-conditional-operator): branchless clamp — the Fixed-vs-GLM perf gate measures this (if/else was 1.5-1.7x slower)1222
return t * t * (T{3} - T{2} * t);1223
});1224
}1226
// ---- Matrix builtins that are not the ordinary product -----------------------------------------1228
/**1229
* The elementwise (Hadamard) product — the GLSL `matrixCompMult`. Named apart from `operator*`1230
* precisely because `*` is the matrix product; this multiplies corresponding entries.1231
* @tparam T the element type. @tparam R the rows. @tparam C the columns.1232
* @param a,b the matrices.1233
* @return the elementwise product.1234
* @complexity O(R·C). @alloc none.1235
* @test Fixarray.MatrixExtras1236
*/1237
template <ndarray::Field T, std::size_t R, std::size_t C>1238
constexpr Mat<T, R, C> matrix_comp_mult(const Mat<T, R, C>& a, const Mat<T, R, C>& b) {1239
return Mat<T, R, C>::from_indices([&](std::size_t i) { return a.data()[i] * b.data()[i]; });1240
}1242
/**1243
* The outer product of a column and a row — the GLSL `outerProduct`: an `R×C` matrix whose1244
* `(i, j)` entry is `c[i] · r[j]`. A rank-one update, the workhorse of a covariance accumulation.1245
* @tparam T the element type. @tparam R the length of @p c (the rows). @tparam C the length of1246
* @p r (the columns).1247
* @param c the column vector. @param r the row vector.1248
* @return the `R×C` outer product.1249
* @complexity O(R·C). @alloc none.1250
* @test Fixarray.MatrixExtras1251
*/1252
template <ndarray::Field T, std::size_t R, std::size_t C>1253
constexpr Mat<T, R, C> outer_product(const Vec<T, R>& c, const Vec<T, C>& r) {1254
// Column-major flat index k addresses row k%R of column k/R, so element k is c[k%R] * r[k/R].1255
return Mat<T, R, C>::from_indices([&](std::size_t k) { return c[k % R] * r[k / R]; });1256
}1258
/**1259
* The inverse transpose of a matrix — `transpose(inverse(m))`, the GLSL `inverseTranspose`. This is1260
* the matrix that carries normals correctly under a non-uniform transform, so lighting stays right.1261
* @tparam T the element type; floating-point.1262
* @tparam N the dimension.1263
* @param m the matrix; must be non-singular.1264
* @return `(m⁻¹)ᵀ`.1265
* @throws std::domain_error when @p m is singular (via @ref inverse).1266
* @complexity O(1) at the fixed sizes. @alloc none.1267
* @test Fixarray.MatrixExtras1268
*/1269
template <ndarray::FloatingPoint T, std::size_t N>1270
constexpr Mat<T, N, N> inverse_transpose(const Mat<T, N, N>& m) {1271
return transpose(inverse(m));1272
}1274
// ---- Named rows and columns: where an enum earns its keep --------------------------------------1275
// GLM indexes a matrix by column (`m[j]`). These free accessors do the same for a @ref Fixed, and1276
// take an @ref ndarray::Subscript — a plain integer OR a scoped `enum class` whose ordinal names the1277
// axis — so a basis vector reads as `column(view, Axis::Forward)` while `Axis` stays a strong type1278
// everywhere else. This is the same door the index operators open, kept open for the free functions.1280
/**1281
* Extract one row of a matrix as a vector.1282
* @tparam T the element type. @tparam R the rows. @tparam C the columns.1283
* @tparam Ix the index type: an integer, or a scoped `enum class` naming the row.1284
* @param m the matrix. @param i the row, `0 <= i < R`.1285
* @return the `C`-vector of that row.1286
* @complexity O(C). @alloc none.1287
* @test Fixarray.NamedRowsAndColumns1288
*/1289
template <ndarray::Field T, std::size_t R, std::size_t C, ::cheatah::ndarray::Subscript Ix>1290
constexpr Vec<T, C> row(const Mat<T, R, C>& m, Ix i) {1291
const auto ri = static_cast<std::size_t>(::cheatah::ndarray::subscript_index(i));1292
Vec<T, C> result;1293
for (std::size_t c = 0; c < C; ++c) { result[c] = m(ri, c); }1294
return result;1295
}1297
/**1298
* Extract one column of a matrix as a vector — a basis vector of the transform. This is the axis a1299
* scoped enum was made to name: `column(view, Axis::Right)`.1300
* @tparam T the element type. @tparam R the rows. @tparam C the columns.1301
* @tparam Ix the index type: an integer, or a scoped `enum class` naming the column.1302
* @param m the matrix. @param j the column, `0 <= j < C`.1303
* @return the `R`-vector of that column.1304
* @complexity O(R). @alloc none.1305
* @test Fixarray.NamedRowsAndColumns1306
*/1307
template <ndarray::Field T, std::size_t R, std::size_t C, ::cheatah::ndarray::Subscript Ix>1308
constexpr Vec<T, R> column(const Mat<T, R, C>& m, Ix j) {1309
const auto cj = static_cast<std::size_t>(::cheatah::ndarray::subscript_index(j));1310
Vec<T, R> result;1311
for (std::size_t r = 0; r < R; ++r) { result[r] = m(r, cj); }1312
return result;1313
}1315
// ---- display ----1316
/**1317
* Render @p v the way an `NDArray` renders — numpy-style nested brackets, each element through the1318
* SHARED scalar formatter (so `i8`/`u8` elements print as NUMBERS, `f32`/`f64` plainly, and a1319
* `complex` as `a+bj`). A vector is `[a, b, c]`; a matrix is `[[…], […]]` in reading `(row, column)`1320
* order — regardless of the column-major storage. This is what `io.print`/`io.str`/`str()` show.1321
* @param v the value to format.1322
* @return the bracketed text.1323
* @complexity O(@ref Fixed::size).1324
* @alloc allocates the result string and a formatting stream per element.1325
* @test Fixarray.ToStringMatchesTheNDArrayRendering1326
*/1327
template <ndarray::Field T, std::size_t... Dims>1328
std::string to_string(const Fixed<T, Dims...>& v) {1329
using F = Fixed<T, Dims...>;1330
std::string out = "[";1331
if constexpr (F::rank == 1) {1332
for (std::size_t i = 0; i < F::size; ++i) {1333
if (i != 0) out += ", ";1334
out += ::cheatah::ndarray::detail::format_scalar(v[i]);1335
}1336
} else {1337
for (std::size_t r = 0; r < F::rows; ++r) {1338
if (r != 0) out += ", ";1339
out += "[";1340
for (std::size_t c = 0; c < F::cols; ++c) {1341
if (c != 0) out += ", ";1342
out += ::cheatah::ndarray::detail::format_scalar(v(r, c));1343
}1344
out += "]";1345
}1346
}1347
return out + "]";1348
}1350
/**1351
* Stream @p v (the nested-bracket @ref to_string form), so a `Fixed` is directly Streamable — a1352
* cheatah `io.print(v)` / `io.str(v)` finds this by ADL, exactly as it does for an `NDArray` or a1353
* primitive.1354
* @param os the stream. @param v the value. @return @p os.1355
* @complexity O(@ref Fixed::size).1356
* @alloc allocates the intermediate string and a formatting stream per element.1357
* @test Fixarray.StreamInsertionUsesTheToStringForm1358
*/1359
template <ndarray::Field T, std::size_t... Dims>1360
std::ostream& operator<<(std::ostream& os, const Fixed<T, Dims...>& v) {1361
return os << to_string(v);1362
}1364
} // namespace cheatah::fixarray1366
// cheatah's value-position subscript `v[i]` / `m[i, j]` lowers to builtins::index(obj, i, ...).1367
// These give it the fixarray meaning: a vector element via operator[], a matrix element via1368
// operator(row, col). An index may be a scoped-enum column label (ndarray::Subscript), matching the1369
// NDArray subscript. (The ndarray overloads live beside these; both are found by the qualified call.)1370
namespace cheatah::builtins {1372
/** Vector element read `v[i]`. @param v the vector. @param i the index (or enum label). @return the element. @complexity O(1). @alloc none. @test Fixarray.BuiltinsIndexLowersSubscripts */1373
template <::cheatah::ndarray::Field T, std::size_t... Dims, ::cheatah::ndarray::Subscript Ix>1374
T index(const ::cheatah::fixarray::Fixed<T, Dims...>& v, Ix i) {1375
return v[i];1376
}1378
/** Matrix element read `m[i, j]`. @param m the matrix. @param i the row. @param j the column (or enum label). @return the element. @complexity O(1). @alloc none. @test Fixarray.BuiltinsIndexLowersSubscripts */1379
template <::cheatah::ndarray::Field T, std::size_t... Dims,1380
::cheatah::ndarray::Subscript I, ::cheatah::ndarray::Subscript J>1381
T index(const ::cheatah::fixarray::Fixed<T, Dims...>& m, I i, J j) {1382
return m(i, j);1383
}1385
/**1386
* Write @p rhs into `v[lo:hi]` — a fixed-size array is FILLED by an assignment, never resized.1387
*1388
* Restricted to `rank == 1`, the same constraint @ref Fixed::operator[] carries: a vector has one1389
* obvious axis to slice, a matrix does not. Bounds follow the list rules (negatives count from the1390
* end, out-of-range clamps, a reversed range is empty). Because the extent is part of the type,1391
* a source of the wrong length is an error rather than a partial write.1392
* @tparam T the element type.1393
* @tparam N the vector's extent.1394
* @tparam R the source range type.1395
* @param v the vector to write into.1396
* @param lo first element (negative counts from the end).1397
* @param hi one past the last, or the "to the end" sentinel.1398
* @param rhs the elements to copy in.1399
* @complexity O(hi - lo).1400
* @alloc none — the destination already owns its storage.1401
* @test Fixarray.SliceAssignCopiesIn1402
* @crtest LangFeatures.FixarraySliceAssignment1403
*/1404
template <::cheatah::ndarray::Field T, std::size_t N, typename R>1405
requires requires(const R& r) { r.begin(); r.end(); }1406
void slice_assign(::cheatah::fixarray::Fixed<T, N>& v, long long lo, long long hi, const R& rhs) {1407
constexpr auto n = static_cast<long long>(N);1408
lo = lo < 0 ? lo + n : lo;1409
if (hi == std::numeric_limits<long long>::max()) {1410
hi = n; // the "to the end" sentinel1411
} else if (hi < 0) {1412
hi += n; // negative counts back from the end1413
}1414
if (lo < 0) lo = 0;1415
if (lo > n) lo = n;1416
if (hi > n) hi = n;1417
if (hi < lo) hi = lo;1418
const auto want = static_cast<std::size_t>(hi - lo);1419
const auto got = static_cast<std::size_t>(std::distance(rhs.begin(), rhs.end()));1420
if (got != want) {1421
throw std::runtime_error(1422
"fixarray: a slice assignment fills a fixed extent — the source has " +1423
std::to_string(got) + " element(s) for " + std::to_string(want) + " slot(s)");1424
}1425
std::copy(rhs.begin(), rhs.end(), v.data() + static_cast<std::size_t>(lo));1426
}1428
} // namespace cheatah::builtins