cheatah
Source BigBrain
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: ndarray
4#pragma once
6/**
7 * @file fixarray.hpp
8 * @brief cheatah `fixarray` — fixed-extent arrays (@ref cheatah::fixarray::Fixed): exactly like an
9 * @ref cheatah::ndarray::NDArray, only faster.
10 *
11 * An @ref cheatah::ndarray::NDArray carries its shape at runtime and its elements on the heap, which
12 * is what makes it general. When the shape is known at compile time and tiny — a 3-D direction, a
13 * 4×4 transform — that generality is the whole cost: a heap allocation, a stride computation and an
14 * 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 are
17 * template parameters, the elements live inline (a `std::array`, so the value is trivially copyable
18 * and sits on the stack or straight inside another struct), the loops have compile-time trip counts
19 * and auto-vectorize, and nothing allocates. **These are the types to reach for in high-performance
20 * applications — a renderer's transforms, a physics solver's contact frames, a filter's small
21 * 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 @ref
24 * cheatah::ndarray::Field, the mathematical index is `(row, column)`, the vocabulary is numpy's
25 * (@ref dot, @ref matmul, @ref transpose, @ref determinant, @ref inverse), and results agree
26 * elementwise. Reach for `NDArray` when the shape is data; reach for `Fixed` when the shape is a
27 * fact about the program.
28 *
29 * **One deliberate difference: a matrix is stored COLUMN-MAJOR**, where `NDArray` is row-major. The
30 * indexing you write is unchanged — `m(row, col)` means what it says, and the constructor still
31 * 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, and
33 * vectorizes, instead of four horizontal dot products that cost a shuffle network; and the buffer is
34 * already in the order graphics APIs (GLSL, SPIR-V, Metal) and GLM expect, so uploading a transform
35 * 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 allocation
40 * mat4f m = mat4f::identity(); // a 4x4, 64 bytes — exactly a push constant
41 * vec3f v = normalize(cross(up, w)); // numpy's vocabulary, glm's speed
42 * ```
43 *
44 * Rank 1 (a vector) and rank 2 (a matrix) are supported; higher ranks are a mechanical extension of
45 * 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's
48 * optimization flags apply. The `linalg` module remains the home of the heavy, shape-generic numerics
49 * on `NDArray` (LU, QR, SVD, eigen); `Fixed` owns the small closed forms where a general factorization
50 * would cost more than the answer.
51 *
52 * **Performance.** Benchmarked against [GLM](https://github.com/g-truc/glm) over the complete overlap
53 * of the two APIs — 160 pairs, every operation, sizes 2/3/4, `float` and `double`, with the outputs
54 * verified identical before either is timed — `Fixed` is **faster than or at parity with GLM on every
55 * one** (19 faster, 141 at parity, none slower; medians over 9 interleaved repetitions, a win counting
56 * only above both 1.15x and 0.25 ns). It wins where structure pays: `mat4f::identity()`
57 * 2.67×, `mat4f * mat4f` 1.72×, `mat4f + mat4f` 1.96×, `inverse(mat4d)` 1.38×. No intrinsics — the code is
58 * 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 <stdexcept>
68#include <string>
69#include <utility>
71#include "ndarray.hpp"
73namespace cheatah::fixarray {
75/// The product of a pack of extents — a fixed array's element count, `1` for an empty pack.
76/// @tparam Dims the extents.
77template <std::size_t... Dims>
78inline constexpr std::size_t extent_product = (std::size_t{1} * ... * Dims);
80/**
81 * Constrains a @ref Fixed to a supported rank: 1 (a vector) or 2 (a matrix). Higher ranks are a
82 * mechanical extension of the same storage, added when a caller needs one — the concept is what
83 * turns "not yet" into a readable compile error instead of a template-instantiation wall.
84 * @tparam Rank the number of extents.
85 */
86template <std::size_t Rank>
87concept SupportedRank = (Rank == 1 || Rank == 2);
89/**
90 * A fixed-extent, inline-stored array — an @ref cheatah::ndarray::NDArray whose shape lives in the
91 * type. Trivially copyable and allocation-free; a matrix is stored column-major (see the file doc).
92 *
93 * @tparam T the element type; any @ref cheatah::ndarray::Field, exactly as `NDArray` accepts.
94 * @tparam Dims the extents. One extent is a vector, two are a matrix (rows, then columns).
95 */
96template <ndarray::Field T, std::size_t... Dims>
97 requires SupportedRank<sizeof...(Dims)> && (((Dims > 0) && ...))
98class Fixed {
99 public:
100 /// The element type.
101 using value_type = T;
103 /// The number of extents: 1 for a vector, 2 for a matrix.
104 static constexpr std::size_t rank = sizeof...(Dims);
105 /// The total number of elements.
106 static constexpr std::size_t size = extent_product<Dims...>;
107 /// The extents, in order (rows, then columns for a matrix).
108 static constexpr std::array<std::size_t, rank> shape{Dims...};
110 /// Rows — the first extent (a vector has one row).
111 static constexpr std::size_t rows = shape[0];
112 /// Columns — the second extent, or 1 for a vector.
113 static constexpr std::size_t cols = rank == 2 ? shape[1] : 1;
115 /// Every element zero — the additive identity, and what a default-constructed value holds.
116 /// @complexity O(size). @alloc none.
117 /// @test Fixarray.DefaultIsZero
118 constexpr Fixed() = default;
120 /**
121 * Construct from exactly @ref size elements, written in READING order: a matrix is given row by
122 * row, the way it appears on paper, regardless of how it is stored. Arguments are converted to
123 * @p T, so a `vec3f` accepts the doubles a cheatah program computes with.
124 * @tparam Args the argument types; each must be convertible to @p T.
125 * @param args the elements in reading order; exactly @ref size of them.
126 * @complexity O(size).
127 * @alloc none.
128 * @test Fixarray.MatrixIndexing
129 * @crtest FixarrayCompileRun.ConstructAndDot
130 */
131 template <class... Args>
132 requires(sizeof...(Args) == size) && (std::convertible_to<Args, T> && ...)
133 explicit constexpr Fixed(Args... args) {
134 const std::array<T, size> reading_order{static_cast<T>(args)...};
135 if constexpr (rank == 1) {
136 data_ = reading_order;
137 } else {
138 for (std::size_t r = 0; r < rows; ++r) {
139 for (std::size_t c = 0; c < cols; ++c) { data_[c * rows + r] = reading_order[r * cols + c]; }
140 }
141 }
142 }
144 /**
145 * The square identity: ones on the diagonal, zeros elsewhere.
146 * @return the identity matrix.
147 * @complexity O(size).
148 * @alloc none.
149 * @test Fixarray.Identity
150 */
151 static constexpr Fixed identity()
152 requires(rank == 2 && rows == cols)
153 {
154 // Built element by element in place. Zeroing the buffer and then poking the diagonal would
155 // store every byte twice; this stores each once, and the compiler folds it to a constant.
156 // In a square column-major buffer the diagonal is exactly the indices divisible by rows + 1.
157 return identity_impl(std::make_index_sequence<size>{});
158 }
160 /**
161 * Every element set to @p value — `filled(0)` is the zero value, `filled(1)` a matrix of ones.
162 * @param value the element to repeat.
163 * @return the filled array.
164 * @complexity O(size).
165 * @alloc none.
166 * @test Fixarray.Filled
167 */
168 static constexpr Fixed filled(T value) {
169 Fixed result;
170 for (std::size_t i = 0; i < size; ++i) { result.data_[i] = value; }
171 return result;
172 }
174 /**
175 * Build an array elementwise: each element `i` of the flat, contiguous buffer is `f(i)`. This is
176 * the allocation-free, single-pass way to write a component-wise operation — no default zeroing
177 * and no separate copy to overwrite, so a call like `abs` or `min` compiles to one vector pass
178 * (`minps`/`maxpd`) rather than two. The index `i` runs over the storage order (column-major for
179 * a matrix), which is exactly what an elementwise operation wants.
180 * @tparam F a callable `T(std::size_t)`.
181 * @param f produces element `i` from its flat index.
182 * @return the array whose element `i` is `f(i)`.
183 * @complexity O(size).
184 * @alloc none.
185 * @test Fixarray.FromIndices
186 */
187 template <class F>
188 static constexpr Fixed from_indices(F&& f) {
189 return from_indices_impl(std::forward<F>(f), std::make_index_sequence<size>{});
190 }
192 /**
193 * Element @p i of a vector.
194 * @param i the index, `0 <= i < size`.
195 * @return a reference to the element.
196 * @complexity O(1).
197 * @alloc none.
198 * @test Fixarray.VectorIndexing
199 */
200 constexpr T& operator[](std::size_t i)
201 requires(rank == 1)
202 {
203 return data_[i];
204 }
206 /// The same, indexed by a scoped `enum class` column label (see @ref ndarray::Subscript): the one
207 /// place an enum is spent as an index, so `v[Axis::Z]` reads column Z while `Axis` stays strong
208 /// everywhere else.
209 /// @tparam Ix the enum index type.
210 /// @param i the element to address, named by an enumerator.
211 /// @return a reference to the element.
212 /// @complexity O(1). @alloc none.
213 /// @test Fixarray.EnumIndexingOnVectorsAndMatrices
214 template <::cheatah::ndarray::Subscript Ix>
215 requires(rank == 1 && std::is_enum_v<Ix>)
216 constexpr T& operator[](Ix i) {
217 return data_[static_cast<std::size_t>(::cheatah::ndarray::subscript_index(i))];
218 }
220 /**
221 * Element @p i of a vector (read-only).
222 * @param i the index, `0 <= i < size`.
223 * @return a const reference to the element.
224 * @complexity O(1).
225 * @alloc none.
226 * @test Fixarray.VectorIndexing
227 */
228 constexpr const T& operator[](std::size_t i) const
229 requires(rank == 1)
230 {
231 return data_[i];
232 }
234 /// Read-only element by a scoped `enum class` column label (see @ref ndarray::Subscript).
235 /// @tparam Ix the enum index type.
236 /// @param i the element to address, named by an enumerator.
237 /// @return a const reference to the element.
238 /// @complexity O(1). @alloc none.
239 /// @test Fixarray.EnumIndexingOnVectorsAndMatrices
240 template <::cheatah::ndarray::Subscript Ix>
241 requires(rank == 1 && std::is_enum_v<Ix>)
242 constexpr const T& operator[](Ix i) const {
243 return data_[static_cast<std::size_t>(::cheatah::ndarray::subscript_index(i))];
244 }
246 /**
247 * Element (@p row, @p col) of a matrix. The index is mathematical; the storage is column-major.
248 * @param row the row, `0 <= row < rows`.
249 * @param col the column, `0 <= col < cols`.
250 * @return a reference to the element.
251 * @complexity O(1).
252 * @alloc none.
253 * @test Fixarray.MatrixIndexing
254 */
255 constexpr T& operator()(std::size_t row, std::size_t col)
256 requires(rank == 2)
257 {
258 return data_[col * rows + row];
259 }
261 /// The same, with either index a scoped `enum class` label (see @ref ndarray::Subscript) — a named
262 /// row or column of a fixed matrix. Mixed integer/enum is allowed; at least one must be an enum, so
263 /// the plain `std::size_t` overload still owns the all-integer call.
264 /// @tparam R the row index type. @tparam C the column index type; at least one is an enum.
265 /// @param row the row to address. @param col the column to address.
266 /// @return a reference to the element.
267 /// @complexity O(1). @alloc none.
268 /// @test Fixarray.EnumIndexingOnVectorsAndMatrices
269 template <::cheatah::ndarray::Subscript R, ::cheatah::ndarray::Subscript C>
270 requires(rank == 2 && (std::is_enum_v<R> || std::is_enum_v<C>))
271 constexpr T& operator()(R row, C col) {
272 return (*this)(static_cast<std::size_t>(::cheatah::ndarray::subscript_index(row)),
273 static_cast<std::size_t>(::cheatah::ndarray::subscript_index(col)));
274 }
276 /**
277 * Element (@p row, @p col) of a matrix, read-only. Mathematical index; column-major storage.
278 * @param row the row, `0 <= row < rows`.
279 * @param col the column, `0 <= col < cols`.
280 * @return a const reference to the element.
281 * @complexity O(1).
282 * @alloc none.
283 * @test Fixarray.MatrixIndexing
284 */
285 constexpr const T& operator()(std::size_t row, std::size_t col) const
286 requires(rank == 2)
287 {
288 return data_[col * rows + row];
289 }
291 /// Read-only (@p row, @p col) with either index a scoped `enum class` label (see @ref
292 /// ndarray::Subscript).
293 /// @tparam R the row index type. @tparam C the column index type; at least one is an enum.
294 /// @param row the row to address. @param col the column to address.
295 /// @return a const reference to the element.
296 /// @complexity O(1). @alloc none.
297 /// @test Fixarray.EnumIndexingOnVectorsAndMatrices
298 template <::cheatah::ndarray::Subscript R, ::cheatah::ndarray::Subscript C>
299 requires(rank == 2 && (std::is_enum_v<R> || std::is_enum_v<C>))
300 constexpr const T& operator()(R row, C col) const {
301 return (*this)(static_cast<std::size_t>(::cheatah::ndarray::subscript_index(row)),
302 static_cast<std::size_t>(::cheatah::ndarray::subscript_index(col)));
303 }
305 /**
306 * A pointer to the elements, contiguous — column-major for a matrix, which is exactly the order a
307 * GPU uniform, a push constant or a BLAS call expects, so an upload is a copy not a transpose.
308 * @return the first element's address.
309 * @complexity O(1).
310 * @alloc none.
311 * @test Fixarray.Data
312 */
313 constexpr T* data() { return data_.data(); }
315 /**
316 * A pointer to the elements, contiguous and column-major for a matrix (read-only).
317 * @return the first element's address.
318 * @complexity O(1).
319 * @alloc none.
320 * @test Fixarray.Data
321 */
322 constexpr const T* data() const { return data_.data(); }
324 /**
325 * Elementwise equality. Exact, as `==` on the elements is exact — floating-point values compare
326 * only if they are bit-for-bit equal.
327 * @param other the array to compare with.
328 * @return true iff every element matches.
329 * @complexity O(size).
330 * @alloc none.
331 * @test Fixarray.Equality
332 */
333 constexpr bool operator==(const Fixed& other) const = default;
335 /**
336 * Add @p other elementwise, in place.
337 * @param other the array to add.
338 * @return a reference to this array.
339 * @complexity O(size).
340 * @alloc none.
341 * @test Fixarray.Arithmetic
342 */
343 constexpr Fixed& operator+=(const Fixed& other) {
344 for (std::size_t i = 0; i < size; ++i) { data_[i] += other.data_[i]; }
345 return *this;
346 }
348 /**
349 * Subtract @p other elementwise, in place.
350 * @param other the array to subtract.
351 * @return a reference to this array.
352 * @complexity O(size).
353 * @alloc none.
354 * @test Fixarray.Arithmetic
355 */
356 constexpr Fixed& operator-=(const Fixed& other) {
357 for (std::size_t i = 0; i < size; ++i) { data_[i] -= other.data_[i]; }
358 return *this;
359 }
361 /**
362 * Scale every element by @p scalar, in place.
363 * @param scalar the factor.
364 * @return a reference to this array.
365 * @complexity O(size).
366 * @alloc none.
367 * @test Fixarray.Arithmetic
368 */
369 constexpr Fixed& operator*=(T scalar) {
370 for (std::size_t i = 0; i < size; ++i) { data_[i] *= scalar; }
371 return *this;
372 }
374 /**
375 * Divide every element by @p scalar, in place.
376 * @param scalar the divisor.
377 * @return a reference to this array.
378 * @complexity O(size).
379 * @alloc none.
380 * @test Fixarray.Arithmetic
381 */
382 constexpr Fixed& operator/=(T scalar) {
383 for (std::size_t i = 0; i < size; ++i) { data_[i] /= scalar; }
384 return *this;
385 }
387 /**
388 * Elementwise sum.
389 * @param a,b the arrays to add.
390 * @return `a + b`.
391 * @complexity O(size).
392 * @alloc none.
393 * @test Fixarray.Arithmetic
394 */
395 friend constexpr Fixed operator+(Fixed a, const Fixed& b) { return a += b; }
397 /**
398 * Elementwise difference.
399 * @param a,b the arrays to subtract.
400 * @return `a - b`.
401 * @complexity O(size).
402 * @alloc none.
403 * @test Fixarray.Arithmetic
404 */
405 friend constexpr Fixed operator-(Fixed a, const Fixed& b) { return a -= b; }
407 /**
408 * Negation.
409 * @param a the array to negate.
410 * @return `-a`.
411 * @complexity O(size).
412 * @alloc none.
413 * @test Fixarray.Arithmetic
414 */
415 friend constexpr Fixed operator-(Fixed a) { return a *= static_cast<T>(-1); }
417 /**
418 * Scale by a scalar.
419 * @param a the array. @param scalar the factor.
420 * @return `a * scalar`.
421 * @complexity O(size).
422 * @alloc none.
423 * @test Fixarray.Arithmetic
424 */
425 friend constexpr Fixed operator*(Fixed a, T scalar) { return a *= scalar; }
427 /**
428 * Scale by a scalar.
429 * @param scalar the factor. @param a the array.
430 * @return `scalar * a`.
431 * @complexity O(size).
432 * @alloc none.
433 * @test Fixarray.Arithmetic
434 */
435 friend constexpr Fixed operator*(T scalar, Fixed a) { return a *= scalar; }
437 /**
438 * Divide by a scalar.
439 * @param a the array. @param scalar the divisor.
440 * @return `a / scalar`.
441 * @complexity O(size).
442 * @alloc none.
443 * @test Fixarray.Arithmetic
444 */
445 friend constexpr Fixed operator/(Fixed a, T scalar) { return a /= scalar; }
447 private:
448 /// Adopt an already-built element buffer, skipping the zero-initialization of the default
449 /// constructor. Private: the buffer's layout (column-major for a matrix) is an implementation
450 /// detail that only the members below may rely on.
451 explicit constexpr Fixed(const std::array<T, size>& elements) : data_(elements) {}
453 /// @ref identity's worker: emits each element exactly once, with no zeroing pass.
454 /// @tparam I the flat indices 0 … size-1.
455 /// @return the identity matrix.
456 template <std::size_t... I>
457 static constexpr Fixed identity_impl(std::index_sequence<I...>)
458 requires(rank == 2 && rows == cols)
459 {
460 return Fixed(std::array<T, size>{(I % (rows + 1) == 0 ? T{1} : T{0})...});
461 }
463 /// @ref from_indices's worker: aggregate-initialises the buffer from `f(0) … f(size-1)`, fully
464 /// unrolled, so there is no loop, no zeroing, and no pointer through which the operands alias.
465 /// @tparam F the element-producing callable.
466 /// @tparam I the flat indices 0 … size-1.
467 /// @param f produces each element from its flat index.
468 /// @return the array of `f(i)`.
469 template <class F, std::size_t... I>
470 static constexpr Fixed from_indices_impl(F&& f, std::index_sequence<I...>) {
471 return Fixed(std::array<T, size>{static_cast<T>(f(I))...});
472 }
474 /// The elements, inline: a vector in order, a matrix column by column. Zero by default.
475 std::array<T, size> data_{};
476};
478/// A fixed-extent vector of @p N elements.
479/// @tparam T the element type. @tparam N the length.
480template <ndarray::Field T, std::size_t N>
481using Vec = Fixed<T, N>;
483/// A fixed-extent matrix of @p R rows and @p C columns — column-major storage, mathematical
484/// `(row, col)` indexing (see the file doc).
485/// @tparam T the element type. @tparam R the rows. @tparam C the columns.
486template <ndarray::Field T, std::size_t R, std::size_t C>
487using Mat = Fixed<T, R, C>;
489/// A 2-D vector of `float`.
490using vec2f = Vec<float, 2>;
491/// A 3-D vector of `float` — a direction, a position, a colour.
492using vec3f = Vec<float, 3>;
493/// A 4-D vector of `float` — a homogeneous point, an RGBA colour.
494using vec4f = Vec<float, 4>;
495/// A 2-D vector of `double`.
496using vec2d = Vec<double, 2>;
497/// A 3-D vector of `double`.
498using vec3d = Vec<double, 3>;
499/// A 4-D vector of `double`.
500using vec4d = Vec<double, 4>;
502/// A 2×2 matrix of `float`.
503using mat2f = Mat<float, 2, 2>;
504/// A 3×3 matrix of `float` — a rotation, or a normal matrix.
505using mat3f = Mat<float, 3, 3>;
506/// A 4×4 matrix of `float` — a transform; exactly the 64 bytes of a push constant.
507using mat4f = Mat<float, 4, 4>;
508/// A 2×2 matrix of `double`.
509using mat2d = Mat<double, 2, 2>;
510/// A 3×3 matrix of `double`.
511using mat3d = Mat<double, 3, 3>;
512/// A 4×4 matrix of `double`.
513using mat4d = Mat<double, 4, 4>;
515namespace detail {
517/**
518 * Sum @p n elements PAIRWISE rather than left to right. A serial `sum += x[i]` chains each add on
519 * the previous one, so the loop runs at the latency of an addition; halving the array instead lets
520 * independent adds issue together, and — the reason numerics people reach for it — the rounding
521 * error grows as O(log n) instead of O(n).
522 * @tparam T the element type. @tparam N the array length.
523 * @param values the products to sum.
524 * @return their sum.
525 * @complexity O(N).
526 * @alloc none.
527 * @test Fixarray.DotAndCross
528 */
529template <ndarray::Field T, std::size_t N>
530constexpr T pairwise_sum(const std::array<T, N>& values) {
531 if constexpr (N == 1) {
532 return values[0];
533 } else if constexpr (N == 2) {
534 return values[0] + values[1];
535 } else if constexpr (N == 3) {
536 return (values[0] + values[1]) + values[2];
537 } else if constexpr (N == 4) {
538 return (values[0] + values[1]) + (values[2] + values[3]);
539 } else {
540 constexpr std::size_t half = N / 2;
541 std::array<T, half> lo{};
542 std::array<T, N - half> hi{};
543 for (std::size_t i = 0; i < half; ++i) { lo[i] = values[i]; }
544 for (std::size_t i = half; i < N; ++i) { hi[i - half] = values[i]; }
545 return pairwise_sum(lo) + pairwise_sum(hi);
546 }
549} // namespace detail
551/**
552 * Inner product of two vectors — Σ aᵢbᵢ, the same quantity @ref dot(const NDArray&, const NDArray&)
553 * computes, without the allocation. Summed pairwise, so it is both faster and more accurate than a
554 * left-to-right accumulation.
555 * @tparam T the element type. @tparam N the length.
556 * @param a,b the vectors.
557 * @return the inner product.
558 * @complexity O(N).
559 * @alloc none.
560 * @test Fixarray.DotAndCross
561 */
562template <ndarray::Field T, std::size_t N>
563constexpr T dot(const Vec<T, N>& a, const Vec<T, N>& b) {
564 // Two regimes, split by what the compiler does with the products. Below 4, an odd width — a
565 // 3-vector of doubles is 24 bytes — spills a products array to the stack, so the terms are
566 // written out and stay in registers. At 4 and above the array is a whole SIMD register (or a
567 // clean multiple), and the loop packs the products into one `mulps`/`mulpd` instead of N scalar
568 // multiplies — which is how this *beats* a scalar dot rather than merely matching it. Both paths
569 // sum pairwise, so the associativity (and the rounding) is identical either way.
570 if constexpr (N == 1) {
571 return a[0] * b[0];
572 } else if constexpr (N == 2) {
573 return a[0] * b[0] + a[1] * b[1];
574 } else if constexpr (N == 3) {
575 return (a[0] * b[0] + a[1] * b[1]) + a[2] * b[2];
576 } else {
577 std::array<T, N> products{};
578 for (std::size_t i = 0; i < N; ++i) { products[i] = a[i] * b[i]; }
579 return detail::pairwise_sum(products);
580 }
583/**
584 * Cross product of two 3-vectors — the vector perpendicular to both, right-handed.
585 * @tparam T the element type.
586 * @param a,b the vectors.
587 * @return `a × b`.
588 * @complexity O(1).
589 * @alloc none.
590 * @test Fixarray.DotAndCross
591 */
592template <ndarray::Field T>
593constexpr Vec<T, 3> cross(const Vec<T, 3>& a, const Vec<T, 3>& b) {
594 return Vec<T, 3>{a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],
595 a[0] * b[1] - a[1] * b[0]};
598/**
599 * The squared Euclidean length of a vector — `dot(v, v)`. Prefer it to @ref norm when only comparing
600 * lengths: it avoids the square root.
601 * @tparam T the element type. @tparam N the length.
602 * @param v the vector.
603 * @return `Σ vᵢ²`.
604 * @complexity O(N).
605 * @alloc none.
606 * @warning for a complex element this is the bilinear `Σ vᵢ²` (matching @ref dot), not
607 * the Hermitian `Σ |vᵢ|²` — it is the squared Euclidean LENGTH only for real
608 * elements.
609 * @test Fixarray.NormAndNormalize
610 */
611template <ndarray::Field T, std::size_t N>
612constexpr T squared_norm(const Vec<T, N>& v) {
613 return dot(v, v);
616/**
617 * The Euclidean length of a vector.
618 * @tparam T the element type; floating-point, since the result is a root.
619 * @tparam N the length.
620 * @param v the vector.
621 * @return `sqrt(Σ vᵢ²)`.
622 * @complexity O(N).
623 * @alloc none.
624 * @test Fixarray.NormAndNormalize
625 */
626template <ndarray::FloatingPoint T, std::size_t N>
627T norm(const Vec<T, N>& v) {
628 return std::sqrt(squared_norm(v));
631/**
632 * The unit vector pointing the same way as @p v.
633 * @tparam T the element type; floating-point.
634 * @tparam N the length.
635 * @param v the vector; must not be the zero vector.
636 * @return `v / norm(v)`.
637 * @throws std::domain_error when @p v has zero length, since it has no direction.
638 * @complexity O(N).
639 * @alloc none.
640 * @test Fixarray.NormAndNormalize
641 */
642template <ndarray::FloatingPoint T, std::size_t N>
643Vec<T, N> normalize(const Vec<T, N>& v) {
644 const T squared = squared_norm(v);
645 if (squared == T{0}) { throw std::domain_error("fixarray::normalize: the zero vector has no direction"); }
646 // One reciprocal, then N multiplies. Dividing each component instead costs N divides, and a
647 // divide is roughly three times the latency of a multiply.
648 const T inverse_length = T{1} / std::sqrt(squared);
649 return v * inverse_length;
652/**
653 * The transpose of a matrix — rows become columns.
654 * @tparam T the element type. @tparam R the rows. @tparam C the columns.
655 * @param m the matrix.
656 * @return the `C×R` transpose.
657 * @complexity O(R·C).
658 * @alloc none.
659 * @test Fixarray.TransposeAndTrace
660 */
661template <ndarray::Field T, std::size_t R, std::size_t C>
662constexpr Mat<T, C, R> transpose(const Mat<T, R, C>& m) {
663 Mat<T, C, R> result;
664 for (std::size_t r = 0; r < R; ++r) {
665 for (std::size_t c = 0; c < C; ++c) { result(c, r) = m(r, c); }
666 }
667 return result;
670/**
671 * The trace of a square matrix — the sum of its diagonal.
672 * @tparam T the element type. @tparam N the dimension.
673 * @param m the matrix.
674 * @return `Σ mᵢᵢ`.
675 * @complexity O(N).
676 * @alloc none.
677 * @test Fixarray.TransposeAndTrace
678 */
679template <ndarray::Field T, std::size_t N>
680constexpr T trace(const Mat<T, N, N>& m) {
681 T sum{};
682 for (std::size_t i = 0; i < N; ++i) { sum += m(i, i); }
683 return sum;
686/**
687 * Matrix product — the same `A·B` @ref matmul(const NDArray&, const NDArray&) computes, with the
688 * shapes checked by the compiler rather than at runtime.
689 * @tparam T the element type. @tparam R the rows of @p a. @tparam K the shared dimension.
690 * @tparam C the columns of @p b.
691 * @param a,b the matrices.
692 * @return the `R×C` product.
693 * @complexity O(R·K·C).
694 * @alloc none.
695 * @test Fixarray.Matmul
696 */
697template <ndarray::Field T, std::size_t R, std::size_t K, std::size_t C>
698constexpr Mat<T, R, C> matmul(const Mat<T, R, K>& a, const Mat<T, K, C>& b) {
699 Mat<T, R, C> result;
700 // Column c of the product is Σₖ (column k of a) · b(k, c): a sum of SCALED COLUMNS. With
701 // column-major storage each of those columns is contiguous, so the inner loop is a plain vertical
702 // multiply-add that vectorizes. The first term SEEDS the column rather than adding to a zeroed
703 // one, which spares a store-then-reload of the accumulator.
704 for (std::size_t c = 0; c < C; ++c) {
705 const T first = b(0, c);
706 for (std::size_t r = 0; r < R; ++r) { result(r, c) = first * a(r, 0); }
707 for (std::size_t k = 1; k < K; ++k) {
708 const T scale = b(k, c);
709 for (std::size_t r = 0; r < R; ++r) { result(r, c) += scale * a(r, k); }
710 }
711 }
712 return result;
715/**
716 * Matrix product, spelled `a * b`.
717 * @tparam T the element type. @tparam R the rows of @p a. @tparam K the shared dimension.
718 * @tparam C the columns of @p b.
719 * @param a,b the matrices.
720 * @return `matmul(a, b)`.
721 * @complexity O(R·K·C).
722 * @alloc none.
723 * @test Fixarray.Matmul
724 */
725template <ndarray::Field T, std::size_t R, std::size_t K, std::size_t C>
726constexpr Mat<T, R, C> operator*(const Mat<T, R, K>& a, const Mat<T, K, C>& b) {
727 return matmul(a, b);
730/**
731 * Transform a vector by a matrix — `A·v`, treating @p v as a column.
732 * @tparam T the element type. @tparam R the rows. @tparam C the columns, and @p v's length.
733 * @param m the matrix. @param v the vector.
734 * @return the `R`-vector `m · v`.
735 * @complexity O(R·C).
736 * @alloc none.
737 * @test Fixarray.Matmul
738 */
739template <ndarray::Field T, std::size_t R, std::size_t C>
740constexpr Vec<T, R> operator*(const Mat<T, R, C>& m, const Vec<T, C>& v) {
741 // `m · v` is Σⱼ (column j of m) · v[j] — a sum of scaled columns, not a stack of row dot
742 // products. Column-major storage makes each column contiguous, so this is a vertical
743 // multiply-add with no horizontal reduction and no shuffles. The first column seeds the result
744 // rather than adding into a zeroed one, which spares a pass over it.
745 Vec<T, R> result;
746 const T first = v[0];
747 for (std::size_t r = 0; r < R; ++r) { result[r] = first * m(r, 0); }
748 for (std::size_t c = 1; c < C; ++c) {
749 const T scale = v[c];
750 for (std::size_t r = 0; r < R; ++r) { result[r] += scale * m(r, c); }
751 }
752 return result;
755/**
756 * The determinant of a 2×2 matrix, in closed form.
757 * @tparam T the element type.
758 * @param m the matrix.
759 * @return `det(m)`.
760 * @complexity O(1).
761 * @alloc none.
762 * @test Fixarray.DeterminantAndInverse
763 */
764template <ndarray::Field T>
765constexpr T determinant(const Mat<T, 2, 2>& m) {
766 return m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0);
769/**
770 * The determinant of a 3×3 matrix, by the rule of Sarrus.
771 * @tparam T the element type.
772 * @param m the matrix.
773 * @return `det(m)`.
774 * @complexity O(1).
775 * @alloc none.
776 * @test Fixarray.DeterminantAndInverse
777 */
778template <ndarray::Field T>
779constexpr T determinant(const Mat<T, 3, 3>& m) {
780 return m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) -
781 m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0)) +
782 m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));
785/**
786 * The determinant of a 4×4 matrix, by cofactor expansion on 2×2 minors — the form a transform
787 * matrix meets, and cheaper than an LU factorization at this size.
788 * @tparam T the element type.
789 * @param m the matrix.
790 * @return `det(m)`.
791 * @complexity O(1).
792 * @alloc none.
793 * @test Fixarray.DeterminantAndInverse
794 */
795template <ndarray::Field T>
796constexpr T determinant(const Mat<T, 4, 4>& m) {
797 const T s0 = m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1);
798 const T s1 = m(0, 0) * m(1, 2) - m(1, 0) * m(0, 2);
799 const T s2 = m(0, 0) * m(1, 3) - m(1, 0) * m(0, 3);
800 const T s3 = m(0, 1) * m(1, 2) - m(1, 1) * m(0, 2);
801 const T s4 = m(0, 1) * m(1, 3) - m(1, 1) * m(0, 3);
802 const T s5 = m(0, 2) * m(1, 3) - m(1, 2) * m(0, 3);
804 const T c5 = m(2, 2) * m(3, 3) - m(3, 2) * m(2, 3);
805 const T c4 = m(2, 1) * m(3, 3) - m(3, 1) * m(2, 3);
806 const T c3 = m(2, 1) * m(3, 2) - m(3, 1) * m(2, 2);
807 const T c2 = m(2, 0) * m(3, 3) - m(3, 0) * m(2, 3);
808 const T c1 = m(2, 0) * m(3, 2) - m(3, 0) * m(2, 2);
809 const T c0 = m(2, 0) * m(3, 1) - m(3, 0) * m(2, 1);
811 return s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0;
814/**
815 * The inverse of a 2×2 matrix, in closed form.
816 * @tparam T the element type; floating-point, since the inverse divides.
817 * @param m the matrix; must be non-singular.
818 * @return `m⁻¹`.
819 * @throws std::domain_error when @p m is singular (zero determinant).
820 * @complexity O(1).
821 * @alloc none.
822 * @test Fixarray.DeterminantAndInverse
823 */
824template <ndarray::FloatingPoint T>
825constexpr Mat<T, 2, 2> inverse(const Mat<T, 2, 2>& m) {
826 const T det = determinant(m);
827 if (det == T{0}) { throw std::domain_error("fixarray::inverse: the matrix is singular"); }
828 const T inv_det = T{1} / det;
829 Mat<T, 2, 2> result;
830 result(0, 0) = m(1, 1) * inv_det;
831 result(0, 1) = -m(0, 1) * inv_det;
832 result(1, 0) = -m(1, 0) * inv_det;
833 result(1, 1) = m(0, 0) * inv_det;
834 return result;
837/**
838 * The inverse of a 3×3 matrix, by its adjugate — the normal matrix a renderer needs.
839 * @tparam T the element type; floating-point.
840 * @param m the matrix; must be non-singular.
841 * @return `m⁻¹`.
842 * @throws std::domain_error when @p m is singular (zero determinant).
843 * @complexity O(1).
844 * @alloc none.
845 * @test Fixarray.DeterminantAndInverse
846 */
847template <ndarray::FloatingPoint T>
848constexpr Mat<T, 3, 3> inverse(const Mat<T, 3, 3>& m) {
849 const T det = determinant(m);
850 if (det == T{0}) { throw std::domain_error("fixarray::inverse: the matrix is singular"); }
851 const T inv_det = T{1} / det;
852 Mat<T, 3, 3> result;
853 result(0, 0) = (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) * inv_det;
854 result(0, 1) = (m(0, 2) * m(2, 1) - m(0, 1) * m(2, 2)) * inv_det;
855 result(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) * inv_det;
856 result(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) * inv_det;
857 result(1, 1) = (m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0)) * inv_det;
858 result(1, 2) = (m(0, 2) * m(1, 0) - m(0, 0) * m(1, 2)) * inv_det;
859 result(2, 0) = (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)) * inv_det;
860 result(2, 1) = (m(0, 1) * m(2, 0) - m(0, 0) * m(2, 1)) * inv_det;
861 result(2, 2) = (m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0)) * inv_det;
862 return result;
865/**
866 * The inverse of a 4×4 matrix, by its adjugate over the 2×2 minors — the transform a camera
867 * inverts every frame.
868 * @tparam T the element type; floating-point.
869 * @param m the matrix; must be non-singular.
870 * @return `m⁻¹`.
871 * @throws std::domain_error when @p m is singular (zero determinant).
872 * @complexity O(1).
873 * @alloc none.
874 * @test Fixarray.DeterminantAndInverse
875 */
876template <ndarray::FloatingPoint T>
877constexpr Mat<T, 4, 4> inverse(const Mat<T, 4, 4>& m) {
878 const T s0 = m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1);
879 const T s1 = m(0, 0) * m(1, 2) - m(1, 0) * m(0, 2);
880 const T s2 = m(0, 0) * m(1, 3) - m(1, 0) * m(0, 3);
881 const T s3 = m(0, 1) * m(1, 2) - m(1, 1) * m(0, 2);
882 const T s4 = m(0, 1) * m(1, 3) - m(1, 1) * m(0, 3);
883 const T s5 = m(0, 2) * m(1, 3) - m(1, 2) * m(0, 3);
885 const T c5 = m(2, 2) * m(3, 3) - m(3, 2) * m(2, 3);
886 const T c4 = m(2, 1) * m(3, 3) - m(3, 1) * m(2, 3);
887 const T c3 = m(2, 1) * m(3, 2) - m(3, 1) * m(2, 2);
888 const T c2 = m(2, 0) * m(3, 3) - m(3, 0) * m(2, 3);
889 const T c1 = m(2, 0) * m(3, 2) - m(3, 0) * m(2, 2);
890 const T c0 = m(2, 0) * m(3, 1) - m(3, 0) * m(2, 1);
892 const T det = s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0;
893 if (det == T{0}) { throw std::domain_error("fixarray::inverse: the matrix is singular"); }
894 const T d = T{1} / det;
896 Mat<T, 4, 4> r;
897 r(0, 0) = (m(1, 1) * c5 - m(1, 2) * c4 + m(1, 3) * c3) * d;
898 r(0, 1) = (-m(0, 1) * c5 + m(0, 2) * c4 - m(0, 3) * c3) * d;
899 r(0, 2) = (m(3, 1) * s5 - m(3, 2) * s4 + m(3, 3) * s3) * d;
900 r(0, 3) = (-m(2, 1) * s5 + m(2, 2) * s4 - m(2, 3) * s3) * d;
902 r(1, 0) = (-m(1, 0) * c5 + m(1, 2) * c2 - m(1, 3) * c1) * d;
903 r(1, 1) = (m(0, 0) * c5 - m(0, 2) * c2 + m(0, 3) * c1) * d;
904 r(1, 2) = (-m(3, 0) * s5 + m(3, 2) * s2 - m(3, 3) * s1) * d;
905 r(1, 3) = (m(2, 0) * s5 - m(2, 2) * s2 + m(2, 3) * s1) * d;
907 r(2, 0) = (m(1, 0) * c4 - m(1, 1) * c2 + m(1, 3) * c0) * d;
908 r(2, 1) = (-m(0, 0) * c4 + m(0, 1) * c2 - m(0, 3) * c0) * d;
909 r(2, 2) = (m(3, 0) * s4 - m(3, 1) * s2 + m(3, 3) * s0) * d;
910 r(2, 3) = (-m(2, 0) * s4 + m(2, 1) * s2 - m(2, 3) * s0) * d;
912 r(3, 0) = (-m(1, 0) * c3 + m(1, 1) * c1 - m(1, 2) * c0) * d;
913 r(3, 1) = (m(0, 0) * c3 - m(0, 1) * c1 + m(0, 2) * c0) * d;
914 r(3, 2) = (-m(3, 0) * s3 + m(3, 1) * s1 - m(3, 2) * s0) * d;
915 r(3, 3) = (m(2, 0) * s3 - m(2, 1) * s1 + m(2, 2) * s0) * d;
916 return r;
919// ---- Geometry: the operations a renderer and a physics solver reach for ------------------------
920// These are the GLSL/GLM geometric builtins, by their standard names, over @ref Fixed vectors: the
921// same mathematics, evaluated in registers with no allocation. They reuse the products above, so a
922// change to @ref dot or the operators reaches them too.
924/**
925 * The Euclidean distance between two points — `norm(a - b)`.
926 * @tparam T the element type; floating-point, since the result is a root.
927 * @tparam N the dimension.
928 * @param a,b the points.
929 * @return `‖a − b‖`.
930 * @complexity O(N). @alloc none.
931 * @test Fixarray.Geometry
932 */
933template <ndarray::FloatingPoint T, std::size_t N>
934T distance(const Vec<T, N>& a, const Vec<T, N>& b) {
935 return norm(a - b);
938/**
939 * The squared distance between two points — `squared_norm(a - b)`. Prefer it to @ref distance when
940 * only comparing distances: it skips the square root.
941 * @tparam T the element type. @tparam N the dimension.
942 * @param a,b the points.
943 * @return `‖a − b‖²`.
944 * @complexity O(N). @alloc none.
945 * @test Fixarray.Geometry
946 */
947template <ndarray::Field T, std::size_t N>
948constexpr T distance_squared(const Vec<T, N>& a, const Vec<T, N>& b) {
949 return squared_norm(a - b);
952/**
953 * Reflect an incident vector about a surface normal — `I − 2 (N·I) N`, the GLSL `reflect`. @p normal
954 * is assumed unit length, as GLSL requires.
955 * @tparam T the element type. @tparam N the dimension.
956 * @param incident the incoming vector.
957 * @param normal the unit surface normal.
958 * @return the reflected vector.
959 * @complexity O(N). @alloc none.
960 * @test Fixarray.Geometry
961 */
962template <ndarray::Field T, std::size_t N>
963constexpr Vec<T, N> reflect(const Vec<T, N>& incident, const Vec<T, N>& normal) {
964 return incident - (T{2} * dot(normal, incident)) * normal;
967/**
968 * Refract an incident vector through a surface — the GLSL `refract`. @p incident and @p normal are
969 * assumed unit length. On total internal reflection (a negative radicand) the result is the zero
970 * vector, exactly as GLSL specifies.
971 * @tparam T the element type; floating-point.
972 * @tparam N the dimension.
973 * @param incident the unit incoming vector.
974 * @param normal the unit surface normal.
975 * @param eta the ratio of indices of refraction (source over destination).
976 * @return the refracted vector, or the zero vector under total internal reflection.
977 * @complexity O(N). @alloc none.
978 * @test Fixarray.Geometry
979 */
980template <ndarray::FloatingPoint T, std::size_t N>
981Vec<T, N> refract(const Vec<T, N>& incident, const Vec<T, N>& normal, T eta) {
982 const T cos_i = dot(normal, incident);
983 const T k = T{1} - eta * eta * (T{1} - cos_i * cos_i);
984 if (k < T{0}) { return Vec<T, N>{}; }
985 return eta * incident - (eta * cos_i + std::sqrt(k)) * normal;
988/**
989 * Orient a normal to face a viewer — the GLSL `faceforward`: return @p n when it already points
990 * against the incident direction, `-n` otherwise. Used to keep a surface normal on the camera's side.
991 * @tparam T the element type. @tparam N the dimension.
992 * @param n the normal to orient.
993 * @param incident the incident vector.
994 * @param reference the reference normal the result is oriented against.
995 * @return @p n or `-n`.
996 * @complexity O(N). @alloc none.
997 * @test Fixarray.Geometry
998 */
999template <ndarray::Numeric T, std::size_t N>
1000constexpr Vec<T, N> faceforward(const Vec<T, N>& n, const Vec<T, N>& incident,
1001 const Vec<T, N>& reference) {
1002 return dot(reference, incident) < T{0} ? n : -n;
1005// ---- Component-wise functions: the GLSL/GLM "common" builtins over a whole array ----------------
1006// Each applies elementwise to every element of a @ref Fixed — a vector or a matrix alike — so a
1007// renderer clamps a colour, a physics step limits a velocity, and a noise field mixes two samples in
1008// the same vocabulary. They read and write the flat buffer, so they are correct whatever the storage
1009// order, and they copy their first argument rather than zero a result and overwrite it.
1011/**
1012 * The absolute value of every element.
1013 * @tparam T the element type; a real number, so `< 0` is meaningful.
1014 * @tparam Dims the extents.
1015 * @param x the array.
1016 * @return `|x|` elementwise.
1017 * @complexity O(size). @alloc none.
1018 * @test Fixarray.CommonUnary
1019 */
1020template <ndarray::Numeric T, std::size_t... Dims>
1021constexpr Fixed<T, Dims...> abs(const Fixed<T, Dims...>& x) {
1022 return Fixed<T, Dims...>::from_indices([&](std::size_t i) {
1023 const T v = x.data()[i];
1024 return v < T{0} ? -v : v;
1025 });
1028/**
1029 * The sign of every element: `-1`, `0`, or `+1`.
1030 * @tparam T the element type; a real number.
1031 * @tparam Dims the extents.
1032 * @param x the array.
1033 * @return the elementwise sign.
1034 * @complexity O(size). @alloc none.
1035 * @test Fixarray.CommonUnary
1036 */
1037template <ndarray::Numeric T, std::size_t... Dims>
1038constexpr Fixed<T, Dims...> sign(const Fixed<T, Dims...>& x) {
1039 return Fixed<T, Dims...>::from_indices([&](std::size_t i) {
1040 const T v = x.data()[i];
1041 return static_cast<T>((T{0} < v) - (v < T{0}));
1042 });
1045/**
1046 * The smaller of each corresponding pair of elements.
1047 * @tparam T the element type; a real number.
1048 * @tparam Dims the extents.
1049 * @param a,b the arrays.
1050 * @return `min(aᵢ, bᵢ)` elementwise.
1051 * @complexity O(size). @alloc none.
1052 * @test Fixarray.MinMaxClamp
1053 */
1054template <ndarray::Numeric T, std::size_t... Dims>
1055constexpr Fixed<T, Dims...> min(const Fixed<T, Dims...>& a, const Fixed<T, Dims...>& b) {
1056 return Fixed<T, Dims...>::from_indices([&](std::size_t i) {
1057 const T ai = a.data()[i];
1058 const T bi = b.data()[i];
1059 return ai < bi ? ai : bi; // a branchless min lowers to minps/minpd
1060 });
1063/**
1064 * Each element floored at the scalar @p s — `min(xᵢ, s)`.
1065 * @tparam T the element type; a real number.
1066 * @tparam Dims the extents.
1067 * @param x the array. @param s the ceiling applied to every element.
1068 * @return `min(xᵢ, s)` elementwise.
1069 * @complexity O(size). @alloc none.
1070 * @test Fixarray.MinMaxClamp
1071 */
1072template <ndarray::Numeric T, std::size_t... Dims>
1073constexpr Fixed<T, Dims...> min(const Fixed<T, Dims...>& x, T s) {
1074 return Fixed<T, Dims...>::from_indices([&](std::size_t i) {
1075 const T v = x.data()[i];
1076 return v < s ? v : s;
1077 });
1080/**
1081 * The larger of each corresponding pair of elements.
1082 * @tparam T the element type; a real number.
1083 * @tparam Dims the extents.
1084 * @param a,b the arrays.
1085 * @return `max(aᵢ, bᵢ)` elementwise.
1086 * @complexity O(size). @alloc none.
1087 * @test Fixarray.MinMaxClamp
1088 */
1089template <ndarray::Numeric T, std::size_t... Dims>
1090constexpr Fixed<T, Dims...> max(const Fixed<T, Dims...>& a, const Fixed<T, Dims...>& b) {
1091 return Fixed<T, Dims...>::from_indices([&](std::size_t i) {
1092 const T ai = a.data()[i];
1093 const T bi = b.data()[i];
1094 return ai < bi ? bi : ai; // branchless max -> maxps/maxpd
1095 });
1098/**
1099 * Each element raised to the scalar @p s — `max(xᵢ, s)`.
1100 * @tparam T the element type; a real number.
1101 * @tparam Dims the extents.
1102 * @param x the array. @param s the floor applied to every element.
1103 * @return `max(xᵢ, s)` elementwise.
1104 * @complexity O(size). @alloc none.
1105 * @test Fixarray.MinMaxClamp
1106 */
1107template <ndarray::Numeric T, std::size_t... Dims>
1108constexpr Fixed<T, Dims...> max(const Fixed<T, Dims...>& x, T s) {
1109 return Fixed<T, Dims...>::from_indices([&](std::size_t i) {
1110 const T v = x.data()[i];
1111 return v < s ? s : v;
1112 });
1115/**
1116 * Constrain every element to `[lo, hi]` — the GLSL `clamp` with scalar bounds, the common case of
1117 * pinning a colour to `[0, 1]`.
1118 * @tparam T the element type; a real number.
1119 * @tparam Dims the extents.
1120 * @param x the array. @param lo the lower bound. @param hi the upper bound.
1121 * @return `min(max(xᵢ, lo), hi)` elementwise.
1122 * @complexity O(size). @alloc none.
1123 * @test Fixarray.MinMaxClamp
1124 */
1125template <ndarray::Numeric T, std::size_t... Dims>
1126constexpr Fixed<T, Dims...> clamp(const Fixed<T, Dims...>& x, T lo, T hi) {
1127 return Fixed<T, Dims...>::from_indices([&](std::size_t i) {
1128 const T v = x.data()[i];
1129 const T low = v < lo ? lo : v;
1130 return hi < low ? hi : low; // min(max(v, lo), hi), branchless
1131 });
1134/**
1135 * Constrain every element between the corresponding bounds — the GLSL `clamp` with per-element
1136 * bounds.
1137 * @tparam T the element type; a real number.
1138 * @tparam Dims the extents.
1139 * @param x the array. @param lo the lower bounds. @param hi the upper bounds.
1140 * @return `min(max(xᵢ, loᵢ), hiᵢ)` elementwise.
1141 * @complexity O(size). @alloc none.
1142 * @test Fixarray.MinMaxClamp
1143 */
1144template <ndarray::Numeric T, std::size_t... Dims>
1145constexpr Fixed<T, Dims...> clamp(const Fixed<T, Dims...>& x, const Fixed<T, Dims...>& lo,
1146 const Fixed<T, Dims...>& hi) {
1147 return Fixed<T, Dims...>::from_indices([&](std::size_t i) {
1148 const T v = x.data()[i];
1149 const T l = lo.data()[i];
1150 const T h = hi.data()[i];
1151 const T low = v < l ? l : v;
1152 return h < low ? h : low;
1153 });
1156/**
1157 * Linear interpolation — the GLSL `mix`: `a (1 − t) + b t`, with a scalar blend @p t (0 gives @p a,
1158 * 1 gives @p b). Composed from the operators, so it inherits their vectorization.
1159 * @tparam T the element type; floating-point.
1160 * @tparam Dims the extents.
1161 * @param a,b the endpoints. @param t the blend factor.
1162 * @return the interpolated array.
1163 * @complexity O(size). @alloc none.
1164 * @test Fixarray.MixStep
1165 */
1166template <ndarray::FloatingPoint T, std::size_t... Dims>
1167constexpr Fixed<T, Dims...> mix(const Fixed<T, Dims...>& a, const Fixed<T, Dims...>& b, T t) {
1168 return a * (T{1} - t) + b * t;
1171/**
1172 * Linear interpolation with a per-element blend — the GLSL `mix` whose factor @p t is an array.
1173 * @tparam T the element type; floating-point.
1174 * @tparam Dims the extents.
1175 * @param a,b the endpoints. @param t the per-element blend factors.
1176 * @return the interpolated array.
1177 * @complexity O(size). @alloc none.
1178 * @test Fixarray.MixStep
1179 */
1180template <ndarray::FloatingPoint T, std::size_t... Dims>
1181constexpr Fixed<T, Dims...> mix(const Fixed<T, Dims...>& a, const Fixed<T, Dims...>& b,
1182 const Fixed<T, Dims...>& t) {
1183 return Fixed<T, Dims...>::from_indices(
1184 [&](std::size_t i) { return a.data()[i] * (T{1} - t.data()[i]) + b.data()[i] * t.data()[i]; });
1187/**
1188 * A step at @p edge — the GLSL `step`: `0` where an element is below @p edge, `1` at or above.
1189 * @tparam T the element type; a real number.
1190 * @tparam Dims the extents.
1191 * @param edge the threshold. @param x the array.
1192 * @return `xᵢ < edge ? 0 : 1` elementwise.
1193 * @complexity O(size). @alloc none.
1194 * @test Fixarray.MixStep
1195 */
1196template <ndarray::Numeric T, std::size_t... Dims>
1197constexpr Fixed<T, Dims...> step(T edge, const Fixed<T, Dims...>& x) {
1198 return Fixed<T, Dims...>::from_indices(
1199 [&](std::size_t i) { return x.data()[i] < edge ? T{0} : T{1}; });
1202/**
1203 * A smooth Hermite transition from 0 to 1 across `[edge0, edge1]` — the GLSL `smoothstep`, with
1204 * everything below @p edge0 giving 0 and everything above @p edge1 giving 1.
1205 * @tparam T the element type; floating-point.
1206 * @tparam Dims the extents.
1207 * @param edge0 the lower edge. @param edge1 the upper edge. @param x the array.
1208 * @return the smoothstepped array.
1209 * @complexity O(size). @alloc none.
1210 * @test Fixarray.MixStep
1211 */
1212template <ndarray::FloatingPoint T, std::size_t... Dims>
1213constexpr Fixed<T, Dims...> smoothstep(T edge0, T edge1, const Fixed<T, Dims...>& x) {
1214 return Fixed<T, Dims...>::from_indices([&](std::size_t i) {
1215 T t = (x.data()[i] - edge0) / (edge1 - edge0);
1216 t = t < T{0} ? T{0} : (T{1} < t ? T{1} : t);
1217 return t * t * (T{3} - T{2} * t);
1218 });
1221// ---- Matrix builtins that are not the ordinary product -----------------------------------------
1223/**
1224 * The elementwise (Hadamard) product — the GLSL `matrixCompMult`. Named apart from `operator*`
1225 * precisely because `*` is the matrix product; this multiplies corresponding entries.
1226 * @tparam T the element type. @tparam R the rows. @tparam C the columns.
1227 * @param a,b the matrices.
1228 * @return the elementwise product.
1229 * @complexity O(R·C). @alloc none.
1230 * @test Fixarray.MatrixExtras
1231 */
1232template <ndarray::Field T, std::size_t R, std::size_t C>
1233constexpr Mat<T, R, C> matrix_comp_mult(const Mat<T, R, C>& a, const Mat<T, R, C>& b) {
1234 return Mat<T, R, C>::from_indices([&](std::size_t i) { return a.data()[i] * b.data()[i]; });
1237/**
1238 * The outer product of a column and a row — the GLSL `outerProduct`: an `R×C` matrix whose
1239 * `(i, j)` entry is `c[i] · r[j]`. A rank-one update, the workhorse of a covariance accumulation.
1240 * @tparam T the element type. @tparam R the length of @p c (the rows). @tparam C the length of
1241 * @p r (the columns).
1242 * @param c the column vector. @param r the row vector.
1243 * @return the `R×C` outer product.
1244 * @complexity O(R·C). @alloc none.
1245 * @test Fixarray.MatrixExtras
1246 */
1247template <ndarray::Field T, std::size_t R, std::size_t C>
1248constexpr Mat<T, R, C> outer_product(const Vec<T, R>& c, const Vec<T, C>& r) {
1249 // Column-major flat index k addresses row k%R of column k/R, so element k is c[k%R] * r[k/R].
1250 return Mat<T, R, C>::from_indices([&](std::size_t k) { return c[k % R] * r[k / R]; });
1253/**
1254 * The inverse transpose of a matrix — `transpose(inverse(m))`, the GLSL `inverseTranspose`. This is
1255 * the matrix that carries normals correctly under a non-uniform transform, so lighting stays right.
1256 * @tparam T the element type; floating-point.
1257 * @tparam N the dimension.
1258 * @param m the matrix; must be non-singular.
1259 * @return `(m⁻¹)ᵀ`.
1260 * @throws std::domain_error when @p m is singular (via @ref inverse).
1261 * @complexity O(1) at the fixed sizes. @alloc none.
1262 * @test Fixarray.MatrixExtras
1263 */
1264template <ndarray::FloatingPoint T, std::size_t N>
1265constexpr Mat<T, N, N> inverse_transpose(const Mat<T, N, N>& m) {
1266 return transpose(inverse(m));
1269// ---- Named rows and columns: where an enum earns its keep --------------------------------------
1270// GLM indexes a matrix by column (`m[j]`). These free accessors do the same for a @ref Fixed, and
1271// take an @ref ndarray::Subscript — a plain integer OR a scoped `enum class` whose ordinal names the
1272// axis — so a basis vector reads as `column(view, Axis::Forward)` while `Axis` stays a strong type
1273// everywhere else. This is the same door the index operators open, kept open for the free functions.
1275/**
1276 * Extract one row of a matrix as a vector.
1277 * @tparam T the element type. @tparam R the rows. @tparam C the columns.
1278 * @tparam Ix the index type: an integer, or a scoped `enum class` naming the row.
1279 * @param m the matrix. @param i the row, `0 <= i < R`.
1280 * @return the `C`-vector of that row.
1281 * @complexity O(C). @alloc none.
1282 * @test Fixarray.NamedRowsAndColumns
1283 */
1284template <ndarray::Field T, std::size_t R, std::size_t C, ::cheatah::ndarray::Subscript Ix>
1285constexpr Vec<T, C> row(const Mat<T, R, C>& m, Ix i) {
1286 const auto ri = static_cast<std::size_t>(::cheatah::ndarray::subscript_index(i));
1287 Vec<T, C> result;
1288 for (std::size_t c = 0; c < C; ++c) { result[c] = m(ri, c); }
1289 return result;
1292/**
1293 * Extract one column of a matrix as a vector — a basis vector of the transform. This is the axis a
1294 * scoped enum was made to name: `column(view, Axis::Right)`.
1295 * @tparam T the element type. @tparam R the rows. @tparam C the columns.
1296 * @tparam Ix the index type: an integer, or a scoped `enum class` naming the column.
1297 * @param m the matrix. @param j the column, `0 <= j < C`.
1298 * @return the `R`-vector of that column.
1299 * @complexity O(R). @alloc none.
1300 * @test Fixarray.NamedRowsAndColumns
1301 */
1302template <ndarray::Field T, std::size_t R, std::size_t C, ::cheatah::ndarray::Subscript Ix>
1303constexpr Vec<T, R> column(const Mat<T, R, C>& m, Ix j) {
1304 const auto cj = static_cast<std::size_t>(::cheatah::ndarray::subscript_index(j));
1305 Vec<T, R> result;
1306 for (std::size_t r = 0; r < R; ++r) { result[r] = m(r, cj); }
1307 return result;
1310// ---- display ----
1311/**
1312 * Render @p v the way an `NDArray` renders — numpy-style nested brackets, each element through the
1313 * SHARED scalar formatter (so `i8`/`u8` elements print as NUMBERS, `f32`/`f64` plainly, and a
1314 * `complex` as `a+bj`). A vector is `[a, b, c]`; a matrix is `[[…], […]]` in reading `(row, column)`
1315 * order — regardless of the column-major storage. This is what `io.print`/`io.str`/`str()` show.
1316 * @param v the value to format.
1317 * @return the bracketed text.
1318 * @complexity O(@ref Fixed::size). @alloc the result string.
1319 * @test Fixarray.ToStringMatchesTheNDArrayRendering
1320 */
1321template <ndarray::Field T, std::size_t... Dims>
1322std::string to_string(const Fixed<T, Dims...>& v) {
1323 using F = Fixed<T, Dims...>;
1324 std::string out = "[";
1325 if constexpr (F::rank == 1) {
1326 for (std::size_t i = 0; i < F::size; ++i) {
1327 if (i != 0) out += ", ";
1328 out += ::cheatah::ndarray::detail::format_scalar(v[i]);
1330 } else {
1331 for (std::size_t r = 0; r < F::rows; ++r) {
1332 if (r != 0) out += ", ";
1333 out += "[";
1334 for (std::size_t c = 0; c < F::cols; ++c) {
1335 if (c != 0) out += ", ";
1336 out += ::cheatah::ndarray::detail::format_scalar(v(r, c));
1338 out += "]";
1341 return out + "]";
1344/**
1345 * Stream @p v (the nested-bracket @ref to_string form), so a `Fixed` is directly Streamable — a
1346 * cheatah `io.print(v)` / `io.str(v)` finds this by ADL, exactly as it does for an `NDArray` or a
1347 * primitive.
1348 * @param os the stream. @param v the value. @return @p os.
1349 * @complexity O(@ref Fixed::size). @alloc the intermediate string.
1350 * @test Fixarray.StreamInsertionUsesTheToStringForm
1351 */
1352template <ndarray::Field T, std::size_t... Dims>
1353std::ostream& operator<<(std::ostream& os, const Fixed<T, Dims...>& v) {
1354 return os << to_string(v);
1357} // namespace cheatah::fixarray
1359// cheatah's value-position subscript `v[i]` / `m[i, j]` lowers to builtins::index(obj, i, ...).
1360// These give it the fixarray meaning: a vector element via operator[], a matrix element via
1361// operator(row, col). An index may be a scoped-enum column label (ndarray::Subscript), matching the
1362// NDArray subscript. (The ndarray overloads live beside these; both are found by the qualified call.)
1363namespace cheatah::builtins {
1365/** 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 */
1366template <::cheatah::ndarray::Field T, std::size_t... Dims, ::cheatah::ndarray::Subscript Ix>
1367T index(const ::cheatah::fixarray::Fixed<T, Dims...>& v, Ix i) {
1368 return v[i];
1371/** 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 */
1372template <::cheatah::ndarray::Field T, std::size_t... Dims,
1373 ::cheatah::ndarray::Subscript I, ::cheatah::ndarray::Subscript J>
1374T index(const ::cheatah::fixarray::Fixed<T, Dims...>& m, I i, J j) {
1375 return m(i, j);
1378} // namespace cheatah::builtins