Source

`stdlib/tests/fixarray_test.cpp`

2// Original work; see ACKNOWLEDGMENTS.md for the open-source ideas we build upon.

3//

4// cheatah::fixarray::Fixed — the fixed-extent arrays. Two things are being proved here:

5//

6// 1. The MATH is right, checked against identities rather than transcribed constants

7// (inverse(m)·m == I, cross(a,b)·a == 0, transpose(transpose(m)) == m). An identity cannot be

8// satisfied by a typo the way a hand-copied expected value can.

9// 2. The COST is right. `Fixed` exists only because NDArray allocates; if it ever stopped being a

10// trivially copyable value of exactly its elements' size, it would have lost its reason to

11// exist. That is a static_assert, not a comment.

12//

13// Both float and double are instantiated: `Fixed` is a template, so an untested instantiation is

14// untested code.

16#include "fixarray.hpp"

18#include <cmath>

19#include <cstdint>

20#include <sstream>

21#include <stdexcept>

22#include <type_traits>

23#include <vector>

25#include <gtest/gtest.h>

27#include "ndarray.hpp"

28#include "routines.hpp"

30namespace fa = cheatah::fixarray;

31namespace la = cheatah::linalg; // linalg routines (inv/det) for the NDArray cross-check below

32namespace nd = cheatah::ndarray;

34namespace {

36/// Elementwise closeness, so a float test and a double test share one predicate.

37template <class M>

38::testing::AssertionResult Close(const M& a, const M& b, double eps = 1e-5) {

39 for (std::size_t i = 0; i < M::size; ++i) {

40 const auto lhs = static_cast<double>(a.data()[i]);

41 const auto rhs = static_cast<double>(b.data()[i]);

42 if (std::fabs(lhs - rhs) > eps) {

43 return ::testing::AssertionFailure()

44 << "element " << i << ": " << lhs << " vs " << rhs;

45 }

46 }

47 return ::testing::AssertionSuccess();

48}

50} // namespace

52// ---- The reason this type exists: no allocation, no padding, no vtable. -------------------------

54TEST(Fixarray, IsAPlainValueOfExactlyItsElements) {

55 static_assert(sizeof(fa::vec2f) == 2 * sizeof(float));

56 static_assert(sizeof(fa::vec3f) == 3 * sizeof(float));

57 static_assert(sizeof(fa::vec4f) == 4 * sizeof(float));

58 static_assert(sizeof(fa::mat3f) == 9 * sizeof(float));

59 static_assert(sizeof(fa::mat4f) == 64); // exactly a 4x4 push constant

60 static_assert(sizeof(fa::mat4d) == 128);

61 static_assert(std::is_trivially_copyable_v<fa::mat4f>);

62 static_assert(std::is_standard_layout_v<fa::mat4f>);

64 // The shape is in the type, so it costs nothing at runtime.

65 static_assert(fa::vec3f::rank == 1);

66 static_assert(fa::vec3f::size == 3);

67 static_assert(fa::vec3f::rows == 3);

68 static_assert(fa::vec3f::cols == 1);

69 static_assert(fa::mat4f::rank == 2);

70 static_assert(fa::mat4f::size == 16);

71 static_assert(fa::mat4f::rows == 4);

72 static_assert(fa::mat4f::cols == 4);

73 static_assert(fa::mat4f::shape[0] == 4 && fa::mat4f::shape[1] == 4);

74 static_assert(std::is_same_v<fa::vec3f::value_type, float>);

75 static_assert(fa::extent_product<2, 3, 4> == 24);

77 // A non-square, non-vector shape is just as ordinary.

78 static_assert(fa::Mat<float, 2, 3>::rows == 2);

79 static_assert(fa::Mat<float, 2, 3>::cols == 3);

80 SUCCEED();

81}

83// ---- Construction, indexing, data() -------------------------------------------------------------

85TEST(Fixarray, DefaultIsZero) {

86 const fa::vec3f v;

87 EXPECT_EQ(v[0], 0.0F);

88 EXPECT_EQ(v[1], 0.0F);

89 EXPECT_EQ(v[2], 0.0F);

90 const fa::mat2d m;

91 EXPECT_EQ(m(0, 0), 0.0);

92 EXPECT_EQ(m(1, 1), 0.0);

93}

95// A fixed-size array is FILLED by a slice assignment — the extent is part of the type, so the

96// values are copied in and nothing is resized.

97TEST(Fixarray, SliceAssignCopiesIn) {

98 using V = fa::Fixed<float, 4>;

99 V v{1.0F, 2.0F, 3.0F, 4.0F};

100 cheatah::builtins::slice_assign(v, 1, 3, std::vector<float>{9.0F, 9.0F});

101 EXPECT_FLOAT_EQ(v[0], 1.0F); // outside the slice: untouched

102 EXPECT_FLOAT_EQ(v[1], 9.0F);

103 EXPECT_FLOAT_EQ(v[2], 9.0F);

104 EXPECT_FLOAT_EQ(v[3], 4.0F);

105 EXPECT_EQ(V::size, 4U); // the extent is compile-time and cannot move

106 // negatives count from the end, exactly as for a list

107 cheatah::builtins::slice_assign(v, -2, -1, std::vector<float>{5.0F});

108 EXPECT_FLOAT_EQ(v[2], 5.0F);

109 // a source of the wrong length is an error, never a partial write

110 EXPECT_THROW(cheatah::builtins::slice_assign(v, 0, 2, std::vector<float>{1.0F}),

111 std::runtime_error);

112 EXPECT_FLOAT_EQ(v[0], 1.0F); // and it really did not write

113}

114

115TEST(Fixarray, VectorIndexing) {

116 fa::vec3f v{1.0F, 2.0F, 3.0F};

117 EXPECT_EQ(v[0], 1.0F);

118 EXPECT_EQ(v[2], 3.0F);

119 v[1] = 9.0F; // non-const

120 EXPECT_EQ(v[1], 9.0F);

121 const fa::vec3f& cv = v;

122 EXPECT_EQ(cv[1], 9.0F); // const

123

124 // Arguments convert: a cheatah program computes in double and stores a float vector.

125 const fa::vec3f from_doubles{1.0, 2.0, 3.0};

126 EXPECT_EQ(from_doubles[2], 3.0F);

127}

128

129TEST(Fixarray, MatrixIndexing) {

130 fa::mat2f m{1.0F, 2.0F, 3.0F, 4.0F}; // row-major

131 EXPECT_EQ(m(0, 0), 1.0F);

132 EXPECT_EQ(m(0, 1), 2.0F);

133 EXPECT_EQ(m(1, 0), 3.0F);

134 EXPECT_EQ(m(1, 1), 4.0F);

135 m(1, 0) = 7.0F; // non-const

136 EXPECT_EQ(m(1, 0), 7.0F);

137 const fa::mat2f& cm = m;

138 EXPECT_EQ(cm(1, 0), 7.0F); // const

139}

140

141TEST(Fixarray, Data) {

142 // The constructor takes elements in READING order...

143 fa::mat2f m{1.0F, 2.0F, 3.0F, 4.0F};

144 EXPECT_EQ(m(0, 0), 1.0F);

145 EXPECT_EQ(m(0, 1), 2.0F);

146 EXPECT_EQ(m(1, 0), 3.0F);

147 EXPECT_EQ(m(1, 1), 4.0F);

148

149 // ...but a matrix is STORED column by column, which is what a GPU uniform, a push constant and

150 // GLM all expect. So the buffer reads 1, 3, 2, 4 — column 0, then column 1.

151 EXPECT_EQ(m.data()[0], 1.0F);

152 EXPECT_EQ(m.data()[1], 3.0F);

153 EXPECT_EQ(m.data()[2], 2.0F);

154 EXPECT_EQ(m.data()[3], 4.0F);

155

156 m.data()[1] = 5.0F; // non-const: element (1, 0), since that is where the buffer says it lives

157 EXPECT_EQ(m(1, 0), 5.0F);

158 const fa::mat2f& cm = m;

159 EXPECT_EQ(cm.data()[3], 4.0F); // const

160

161 // A vector has one order and no ambiguity.

162 const fa::vec3f v{7.0F, 8.0F, 9.0F};

163 EXPECT_EQ(v.data()[1], 8.0F);

164}

165

166TEST(Fixarray, Identity) {

167 constexpr fa::mat3f compile_time = fa::mat3f::identity(); // usable at compile time

168 static_assert(compile_time(0, 0) == 1.0F);

169 static_assert(compile_time(0, 1) == 0.0F);

170

171 // ...and at run time. A constexpr function nobody executes is a function nobody proved runs.

172 fa::mat3f runtime = fa::mat3f::identity();

173 for (std::size_t r = 0; r < 3; ++r) {

174 for (std::size_t c = 0; c < 3; ++c) {

175 EXPECT_EQ(runtime(r, c), r == c ? 1.0F : 0.0F);

176 }

177 }

178 runtime(2, 2) = 5.0F; // the non-const matrix accessor on this instantiation

179 EXPECT_EQ(runtime(2, 2), 5.0F);

180 EXPECT_EQ(fa::mat4d::identity()(3, 3), 1.0);

181 EXPECT_EQ(fa::mat2d::identity()(0, 0), 1.0);

182}

183

184TEST(Fixarray, Filled) {

185 const fa::mat2f threes = fa::mat2f::filled(3.0F);

186 EXPECT_EQ(threes(0, 0), 3.0F);

187 EXPECT_EQ(threes(1, 1), 3.0F);

188 EXPECT_EQ(fa::vec4d::filled(-1.0)[3], -1.0);

189}

190

191TEST(Fixarray, Equality) {

192 const fa::vec3f a{1.0F, 2.0F, 3.0F};

193 const fa::vec3f b{1.0F, 2.0F, 3.0F};

194 const fa::vec3f c{1.0F, 2.0F, 4.0F};

195 EXPECT_TRUE(a == b);

196 EXPECT_FALSE(a == c);

197 EXPECT_TRUE(a != c);

198 EXPECT_FALSE(a != b);

199}

200

201// ---- Arithmetic ---------------------------------------------------------------------------------

202

203TEST(Fixarray, Arithmetic) {

204 const fa::vec3f a{1.0F, 2.0F, 3.0F};

205 const fa::vec3f b{4.0F, 5.0F, 6.0F};

206

207 EXPECT_TRUE(a + b == (fa::vec3f{5.0F, 7.0F, 9.0F}));

208 EXPECT_TRUE(b - a == (fa::vec3f{3.0F, 3.0F, 3.0F}));

209 EXPECT_TRUE(-a == (fa::vec3f{-1.0F, -2.0F, -3.0F}));

210 EXPECT_TRUE(a * 2.0F == (2.0F * a)); // scalar multiply, both orders

211 EXPECT_TRUE((a * 2.0F) / 2.0F == a); // and its inverse

212 EXPECT_TRUE(a + (-a) == fa::vec3f{}); // additive inverse

213

214 fa::vec3f m = a;

215 m += b;

216 EXPECT_TRUE(m == a + b);

217 m -= b;

218 EXPECT_TRUE(m == a);

219 m *= 3.0F;

220 EXPECT_TRUE(m == a * 3.0F);

221 m /= 3.0F;

222 EXPECT_TRUE(m == a);

223

224 // Doubles behave the same.

225 fa::mat2d dm{1.0, 2.0, 3.0, 4.0};

226 dm += fa::mat2d::filled(1.0);

227 EXPECT_EQ(dm(0, 0), 2.0);

228 dm -= fa::mat2d::filled(1.0);

229 EXPECT_EQ(dm(0, 0), 1.0);

230 dm *= 2.0;

231 EXPECT_EQ(dm(1, 1), 8.0);

232 dm /= 2.0;

233 EXPECT_EQ(dm(1, 1), 4.0);

234 EXPECT_EQ((-dm)(1, 1), -4.0);

235 EXPECT_EQ((dm + dm)(0, 0), 2.0);

236 EXPECT_EQ((dm - dm)(0, 0), 0.0);

237 EXPECT_EQ((2.0 * dm)(0, 0), 2.0);

238 EXPECT_EQ((dm / 2.0)(1, 1), 2.0);

239

240 // Every alias is a real instantiation; exercise the smaller ones so none is merely declared.

241 fa::vec2d small{2.0, 4.0};

242 small /= 2.0;

243 EXPECT_EQ(small[1], 2.0);

244 EXPECT_EQ((small / 2.0)[0], 0.5);

245 EXPECT_EQ((fa::vec2f{1.0F, 2.0F} + fa::vec2f{1.0F, 1.0F})[1], 3.0F);

246 EXPECT_EQ((fa::vec4d::filled(2.0) * 0.5)[0], 1.0);

247}

248

249// ---- Vector products ----------------------------------------------------------------------------

250

251TEST(Fixarray, DotAndCross) {

252 constexpr fa::vec3f a{1.0F, 2.0F, 3.0F};

253 constexpr fa::vec3f b{4.0F, 5.0F, 6.0F};

254 static_assert(fa::dot(a, b) == 32.0F); // compile-time

255 EXPECT_EQ(fa::dot(a, b), 32.0F);

256

257 constexpr fa::vec3f c = fa::cross(a, b);

258 static_assert(c[0] == -3.0F && c[1] == 6.0F && c[2] == -3.0F);

259

260 // The identity that defines a cross product: perpendicular to both operands.

261 EXPECT_EQ(fa::dot(c, a), 0.0F);

262 EXPECT_EQ(fa::dot(c, b), 0.0F);

263 // ...and anticommutative.

264 EXPECT_TRUE(fa::cross(b, a) == -c);

265

266 // Right-handed: x cross y == z.

267 const fa::vec3d x{1.0, 0.0, 0.0};

268 const fa::vec3d y{0.0, 1.0, 0.0};

269 EXPECT_TRUE(fa::cross(x, y) == (fa::vec3d{0.0, 0.0, 1.0}));

270 EXPECT_EQ(fa::dot(fa::vec4f{1.0F, 1.0F, 1.0F, 1.0F}, fa::vec4f{1.0F, 2.0F, 3.0F, 4.0F}), 10.0F);

271}

272

273TEST(Fixarray, NormAndNormalize) {

274 const fa::vec3f v{3.0F, 4.0F, 0.0F};

275 EXPECT_EQ(fa::squared_norm(v), 25.0F);

276 EXPECT_FLOAT_EQ(fa::norm(v), 5.0F);

277

278 const fa::vec3f unit = fa::normalize(v);

279 EXPECT_FLOAT_EQ(fa::norm(unit), 1.0F);

280 EXPECT_FLOAT_EQ(unit[0], 0.6F);

281 EXPECT_FLOAT_EQ(unit[1], 0.8F);

282

283 EXPECT_DOUBLE_EQ(fa::norm(fa::vec2d{0.0, 2.0}), 2.0);

284

285 // Every instantiation must normalize, not merely refuse to: `normalize` is one reciprocal and a

286 // multiply, and a size that only ever saw the throw path is a size nobody proved works.

287 EXPECT_DOUBLE_EQ(fa::norm(fa::normalize(fa::vec2d{3.0, 4.0})), 1.0);

288 EXPECT_FLOAT_EQ(fa::norm(fa::normalize(fa::vec2f{0.0F, 2.0F})), 1.0F);

289 EXPECT_DOUBLE_EQ(fa::norm(fa::normalize(fa::vec4d{1.0, 1.0, 1.0, 1.0})), 1.0);

290 EXPECT_FLOAT_EQ(fa::normalize(fa::vec2f{0.0F, 2.0F})[1], 1.0F);

291

292 // The zero vector has no direction; saying so beats returning NaNs.

293 EXPECT_THROW((void)fa::normalize(fa::vec3f{}), std::domain_error);

294 EXPECT_THROW((void)fa::normalize(fa::vec2d{}), std::domain_error);

295 EXPECT_THROW((void)fa::normalize(fa::vec4d{}), std::domain_error);

296}

297

298// ---- Matrix products, transpose, trace ----------------------------------------------------------

299

300TEST(Fixarray, Matmul) {

301 const fa::mat2f a{1.0F, 2.0F, 3.0F, 4.0F};

302 const fa::mat2f b{5.0F, 6.0F, 7.0F, 8.0F};

303 const fa::mat2f ab = fa::matmul(a, b);

304 EXPECT_EQ(ab(0, 0), 19.0F);

305 EXPECT_EQ(ab(0, 1), 22.0F);

306 EXPECT_EQ(ab(1, 0), 43.0F);

307 EXPECT_EQ(ab(1, 1), 50.0F);

308 EXPECT_TRUE(a * b == ab); // the operator spelling

309

310 // Identity is the multiplicative identity, and matmul is associative.

311 EXPECT_TRUE(a * fa::mat2f::identity() == a);

312 EXPECT_TRUE(fa::mat2f::identity() * a == a);

313 const fa::mat2f c{2.0F, 0.0F, 1.0F, 3.0F};

314 EXPECT_TRUE(Close((a * b) * c, a * (b * c)));

315

316 // Non-square shapes chain: (2x3)(3x2) -> 2x2.

317 const fa::Mat<double, 2, 3> wide{1.0, 2.0, 3.0, 4.0, 5.0, 6.0};

318 const fa::Mat<double, 3, 2> tall{7.0, 8.0, 9.0, 10.0, 11.0, 12.0};

319 const fa::Mat<double, 2, 2> product = wide * tall;

320 EXPECT_DOUBLE_EQ(product(0, 0), 58.0);

321 EXPECT_DOUBLE_EQ(product(1, 1), 154.0);

322

323 // Matrix times vector.

324 const fa::vec2f v = a * fa::vec2f{1.0F, 1.0F};

325 EXPECT_EQ(v[0], 3.0F);

326 EXPECT_EQ(v[1], 7.0F);

327 const fa::vec3d w = fa::mat3d::identity() * fa::vec3d{1.0, 2.0, 3.0};

328 EXPECT_TRUE(w == (fa::vec3d{1.0, 2.0, 3.0}));

329 const fa::Vec<double, 2> rect = wide * fa::vec3d{1.0, 1.0, 1.0};

330 EXPECT_DOUBLE_EQ(rect[0], 6.0);

331 EXPECT_DOUBLE_EQ(rect[1], 15.0);

332}

333

334TEST(Fixarray, TransposeAndTrace) {

335 const fa::mat2f m{1.0F, 2.0F, 3.0F, 4.0F};

336 const fa::mat2f t = fa::transpose(m);

337 EXPECT_EQ(t(0, 1), 3.0F);

338 EXPECT_EQ(t(1, 0), 2.0F);

339 EXPECT_TRUE(fa::transpose(t) == m); // an involution

340 EXPECT_EQ(fa::trace(m), 5.0F);

341 EXPECT_EQ(fa::trace(fa::mat4d::identity()), 4.0);

342

343 // A non-square transpose swaps the shape.

344 const fa::Mat<float, 2, 3> wide{1.0F, 2.0F, 3.0F, 4.0F, 5.0F, 6.0F};

345 const fa::Mat<float, 3, 2> narrow = fa::transpose(wide);

346 static_assert(decltype(narrow)::rows == 3 && decltype(narrow)::cols == 2);

347 EXPECT_EQ(narrow(2, 0), 3.0F);

348 EXPECT_EQ(narrow(0, 1), 4.0F);

349}

350

351// ---- Determinant and inverse --------------------------------------------------------------------

352

353TEST(Fixarray, DeterminantAndInverse) {

354 // 2x2

355 const fa::mat2f m2{4.0F, 7.0F, 2.0F, 6.0F};

356 EXPECT_FLOAT_EQ(fa::determinant(m2), 10.0F);

357 EXPECT_TRUE(Close(fa::inverse(m2) * m2, fa::mat2f::identity()));

358 EXPECT_TRUE(Close(m2 * fa::inverse(m2), fa::mat2f::identity()));

359

360 // 3x3

361 const fa::mat3d m3{2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0};

362 EXPECT_DOUBLE_EQ(fa::determinant(m3), 4.0);

363 EXPECT_TRUE(Close(fa::inverse(m3) * m3, fa::mat3d::identity(), 1e-12));

364

365 // 4x4

366 const fa::mat4d m4{1.0, 2.0, 0.0, 1.0, 0.0, 1.0, 3.0, 0.0,

367 2.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0};

368 EXPECT_TRUE(Close(fa::inverse(m4) * m4, fa::mat4d::identity(), 1e-12));

369 EXPECT_TRUE(Close(m4 * fa::inverse(m4), fa::mat4d::identity(), 1e-12));

370 EXPECT_TRUE(Close(fa::inverse(fa::inverse(m4)), m4, 1e-10)); // an involution

371

372 // det(I) == 1 at every supported size, and det(AB) == det(A)det(B).

373 EXPECT_FLOAT_EQ(fa::determinant(fa::mat2f::identity()), 1.0F);

374 EXPECT_DOUBLE_EQ(fa::determinant(fa::mat3d::identity()), 1.0);

375 EXPECT_DOUBLE_EQ(fa::determinant(fa::mat4d::identity()), 1.0);

376 const fa::mat3d other{1.0, 2.0, 3.0, 0.0, 1.0, 4.0, 5.0, 6.0, 0.0};

377 EXPECT_NEAR(fa::determinant(m3 * other), fa::determinant(m3) * fa::determinant(other), 1e-9);

378 EXPECT_FLOAT_EQ(fa::determinant(fa::mat4f{1.0F, 2.0F, 0.0F, 1.0F, 0.0F, 1.0F, 3.0F, 0.0F,

379 2.0F, 0.0F, 1.0F, 1.0F, 1.0F, 1.0F, 1.0F, 2.0F}),

380 static_cast<float>(fa::determinant(m4)));

381

382 // A singular matrix has no inverse, and says so rather than returning infinities.

383 EXPECT_EQ(fa::determinant(fa::mat2f{}), 0.0F);

384 EXPECT_THROW((void)fa::inverse(fa::mat2f{}), std::domain_error);

385 EXPECT_THROW((void)fa::inverse(fa::mat3d{}), std::domain_error);

386 EXPECT_THROW((void)fa::inverse(fa::mat4d{}), std::domain_error);

387

388 // Rank-deficient, not merely all-zero: two identical rows.

389 EXPECT_THROW((void)fa::inverse(fa::mat3d{1.0, 2.0, 3.0, 1.0, 2.0, 3.0, 4.0, 5.0, 7.0}),

390 std::domain_error);

391}

392

393// ---- The type is not secretly limited to the graphics sizes -----------------------------------

394

395TEST(Fixarray, WorksBeyondTheAliasedSizes) {

396 // The aliases stop at 4 because that is where graphics stops; the TYPE does not. This also

397 // exercises `dot`'s general recursive pairwise sum, which the 2/3/4 cases short-circuit past.

398 fa::Vec<double, 8> a;

399 fa::Vec<double, 8> b;

400 for (std::size_t i = 0; i < 8; ++i) {

401 a[i] = static_cast<double>(i + 1); // 1..8

402 b[i] = 1.0;

403 }

404 EXPECT_DOUBLE_EQ(fa::dot(a, b), 36.0); // 1+2+...+8

405 EXPECT_DOUBLE_EQ(fa::squared_norm(b), 8.0); // eight ones

406 EXPECT_DOUBLE_EQ(fa::norm(b), std::sqrt(8.0));

407 EXPECT_DOUBLE_EQ(fa::norm(fa::normalize(a)), 1.0);

408

409 // An odd length exercises the uneven split of the recursion (5 = 2 + 3).

410 fa::Vec<float, 5> odd{1.0F, 2.0F, 3.0F, 4.0F, 5.0F};

411 EXPECT_FLOAT_EQ(fa::dot(odd, odd), 55.0F); // 1+4+9+16+25

412

413 // And a bigger matrix still multiplies, transposes and transforms.

414 const fa::Mat<double, 5, 5> identity5 = fa::Mat<double, 5, 5>::identity();

415 EXPECT_DOUBLE_EQ(fa::trace(identity5), 5.0);

416 EXPECT_TRUE(fa::transpose(identity5) == identity5);

417 const fa::Vec<double, 5> five{1.0, 2.0, 3.0, 4.0, 5.0};

418 const fa::Vec<double, 5> through = identity5 * five;

419 EXPECT_TRUE(through == five);

420 EXPECT_TRUE(fa::matmul(identity5, identity5) == identity5);

421}

422

423// ---- The same answers as NDArray, which is the promise the name makes -------------------------

424

425TEST(Fixarray, AgreesWithTheDynamicNDArray) {

426 // `Fixed` claims to be "NDArray, only faster". The claim is only worth making if the answers

427 // agree; a 3x3 inverse and determinant are where that is easiest to check.

428 const std::vector<double> values{2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0};

429

430 const fa::mat3d fixed{values[0], values[1], values[2], values[3], values[4],

431 values[5], values[6], values[7], values[8]};

432 const fa::mat3d fixed_inv = fa::inverse(fixed);

433

434 const nd::NDArray dynamic = nd::reshape(nd::array(values), {3, 3});

435 const nd::NDArray dynamic_inv = la::inv(dynamic);

436

437 for (long long r = 0; r < 3; ++r) {

438 for (long long c = 0; c < 3; ++c) {

439 EXPECT_NEAR(fixed_inv(static_cast<std::size_t>(r), static_cast<std::size_t>(c)),

440 nd::get(dynamic_inv, {r, c}), 1e-12);

441 }

442 }

443 EXPECT_NEAR(fa::determinant(fixed), la::det(dynamic), 1e-12);

444}

445

446// ---- The GLSL/GLM surface: geometry ------------------------------------------------------------

447

448TEST(Fixarray, Geometry) {

449 const fa::vec3f a{1.0F, 2.0F, 3.0F};

450 const fa::vec3f b{4.0F, 6.0F, 8.0F};

451 EXPECT_FLOAT_EQ(fa::distance(a, b), std::sqrt(9.0F + 16.0F + 25.0F));

452 EXPECT_FLOAT_EQ(fa::distance_squared(a, b), 50.0F);

453 EXPECT_DOUBLE_EQ(fa::distance(fa::vec2d{0.0, 0.0}, fa::vec2d{3.0, 4.0}), 5.0);

454

455 // reflect off the floor (normal +y): the y-component flips, x and z survive.

456 const fa::vec3f down{1.0F, -1.0F, 0.0F};

457 const fa::vec3f up{0.0F, 1.0F, 0.0F};

458 EXPECT_TRUE(fa::reflect(down, up) == (fa::vec3f{1.0F, 1.0F, 0.0F}));

459 // A vector reflected twice about the same normal returns to itself (dot with unit normal).

460 EXPECT_TRUE(fa::reflect(fa::reflect(down, up), up) == down);

461

462 // refract with equal indices (eta = 1) does not bend, so a unit vector stays unit.

463 const fa::vec3f incident = fa::normalize(fa::vec3f{1.0F, -1.0F, 0.0F});

464 EXPECT_FLOAT_EQ(fa::norm(fa::refract(incident, up, 1.0F)), 1.0F);

465 // Total internal reflection returns the zero vector.

466 const fa::vec2d grazing = fa::normalize(fa::vec2d{1.0, -0.01});

467 EXPECT_TRUE(fa::refract(grazing, fa::vec2d{0.0, 1.0}, 5.0) == fa::vec2d{});

468

469 // faceforward keeps a normal on the incident's side. dot(nref, I) < 0 -> return n unchanged.

470 EXPECT_TRUE(fa::faceforward(up, down, up) == up);

471 const fa::vec3f away{0.0F, 1.0F, 0.0F};

472 EXPECT_TRUE(fa::faceforward(up, fa::vec3f{0.0F, 1.0F, 0.0F}, away) == (fa::vec3f{0.0F, -1.0F, 0.0F}));

473}

474

475// ---- The GLSL/GLM surface: component-wise common builtins ---------------------------------------

476

477TEST(Fixarray, CommonUnary) {

478 EXPECT_TRUE(fa::abs(fa::vec4f{-1.0F, 2.0F, -3.0F, 0.0F}) == (fa::vec4f{1.0F, 2.0F, 3.0F, 0.0F}));

479 EXPECT_TRUE(fa::sign(fa::vec3f{-2.0F, 0.0F, 5.0F}) == (fa::vec3f{-1.0F, 0.0F, 1.0F}));

480 // Works on a matrix too — it is elementwise over the whole array.

481 EXPECT_TRUE(fa::abs(fa::mat2d{-1.0, 2.0, -3.0, 4.0}) == (fa::mat2d{1.0, 2.0, 3.0, 4.0}));

482 EXPECT_TRUE(fa::sign(fa::mat2f{-4.0F, 0.0F, 8.0F, -1.0F}) == (fa::mat2f{-1.0F, 0.0F, 1.0F, -1.0F}));

483 EXPECT_TRUE(fa::abs(fa::vec2d{-1.5, -2.5}) == (fa::vec2d{1.5, 2.5}));

484}

485

486TEST(Fixarray, MinMaxClamp) {

487 const fa::vec3f a{1.0F, 5.0F, 3.0F};

488 const fa::vec3f b{4.0F, 2.0F, 6.0F};

489 EXPECT_TRUE(fa::min(a, b) == (fa::vec3f{1.0F, 2.0F, 3.0F}));

490 EXPECT_TRUE(fa::max(a, b) == (fa::vec3f{4.0F, 5.0F, 6.0F}));

491 EXPECT_TRUE(fa::min(fa::vec3f{1.0F, 5.0F, 9.0F}, 4.0F) == (fa::vec3f{1.0F, 4.0F, 4.0F}));

492 EXPECT_TRUE(fa::max(fa::vec3f{1.0F, 5.0F, 9.0F}, 4.0F) == (fa::vec3f{4.0F, 5.0F, 9.0F}));

493

494 // scalar-bound clamp — pinning a colour to [0, 1]

495 EXPECT_TRUE(fa::clamp(fa::vec4f{-1.0F, 0.5F, 2.0F, 0.0F}, 0.0F, 1.0F) ==

496 (fa::vec4f{0.0F, 0.5F, 1.0F, 0.0F}));

497 // per-element bounds

498 EXPECT_TRUE(fa::clamp(fa::vec3d{5.0, -5.0, 0.5}, fa::vec3d{0.0, 0.0, 0.0},

499 fa::vec3d{1.0, 1.0, 1.0}) == (fa::vec3d{1.0, 0.0, 0.5}));

500 // matrices too (double, to exercise that instantiation)

501 EXPECT_TRUE(fa::min(fa::mat2d{1.0, 4.0, 3.0, 2.0}, fa::mat2d{2.0, 2.0, 2.0, 2.0}) ==

502 (fa::mat2d{1.0, 2.0, 2.0, 2.0}));

503 EXPECT_TRUE(fa::max(fa::mat2d::filled(1.0), 3.0) == fa::mat2d::filled(3.0));

504}

505

506TEST(Fixarray, MixStep) {

507 // mix with a scalar factor is a lerp

508 EXPECT_TRUE(fa::mix(fa::vec3f{0.0F, 0.0F, 0.0F}, fa::vec3f{2.0F, 4.0F, 6.0F}, 0.5F) ==

509 (fa::vec3f{1.0F, 2.0F, 3.0F}));

510 EXPECT_TRUE(fa::mix(fa::vec2d{1.0, 1.0}, fa::vec2d{3.0, 5.0}, 0.0) == (fa::vec2d{1.0, 1.0}));

511 EXPECT_TRUE(fa::mix(fa::vec2d{1.0, 1.0}, fa::vec2d{3.0, 5.0}, 1.0) == (fa::vec2d{3.0, 5.0}));

512 // per-element factor

513 EXPECT_TRUE(fa::mix(fa::vec3f{0.0F, 0.0F, 0.0F}, fa::vec3f{10.0F, 10.0F, 10.0F},

514 fa::vec3f{0.0F, 0.5F, 1.0F}) == (fa::vec3f{0.0F, 5.0F, 10.0F}));

515

516 // step: below the edge is 0, at or above is 1

517 EXPECT_TRUE(fa::step(2.0F, fa::vec3f{1.0F, 2.0F, 3.0F}) == (fa::vec3f{0.0F, 1.0F, 1.0F}));

518 EXPECT_TRUE(fa::step(0.0, fa::vec2d{-1.0, 1.0}) == (fa::vec2d{0.0, 1.0}));

519

520 // smoothstep: clamped at the edges, 0.5 at the midpoint, Hermite in between

521 const fa::vec4f s = fa::smoothstep(0.0F, 1.0F, fa::vec4f{-1.0F, 0.0F, 0.5F, 2.0F});

522 EXPECT_FLOAT_EQ(s[0], 0.0F);

523 EXPECT_FLOAT_EQ(s[1], 0.0F);

524 EXPECT_FLOAT_EQ(s[2], 0.5F);

525 EXPECT_FLOAT_EQ(s[3], 1.0F);

526 // monotone and within [0,1]

527 const fa::vec2d q = fa::smoothstep(0.0, 10.0, fa::vec2d{2.5, 7.5});

528 EXPECT_GT(q[1], q[0]);

529 EXPECT_GE(q[0], 0.0);

530 EXPECT_LE(q[1], 1.0);

531}

532

533// ---- The GLSL/GLM surface: matrix builtins -----------------------------------------------------

534

535TEST(Fixarray, MatrixExtras) {

536 // Hadamard product multiplies corresponding entries (NOT the matrix product).

537 const fa::mat2f m{1.0F, 2.0F, 3.0F, 4.0F};

538 const fa::mat2f k{2.0F, 0.0F, 0.0F, 2.0F};

539 EXPECT_TRUE(fa::matrix_comp_mult(m, k) == (fa::mat2f{2.0F, 0.0F, 0.0F, 8.0F}));

540

541 // outer product: (i, j) = c[i] * r[j]

542 const fa::Mat<float, 2, 3> op = fa::outer_product(fa::vec2f{1.0F, 2.0F}, fa::vec3f{3.0F, 4.0F, 5.0F});

543 EXPECT_FLOAT_EQ(op(0, 0), 3.0F);

544 EXPECT_FLOAT_EQ(op(0, 2), 5.0F);

545 EXPECT_FLOAT_EQ(op(1, 1), 8.0F);

546 // outer_product(c, r) == c as a column times r as a row, so its rank is one: rows are multiples.

547 EXPECT_FLOAT_EQ(op(1, 0) / op(0, 0), 2.0F);

548

549 // inverse_transpose carries normals: for an orthonormal (rotation) matrix it equals the matrix

550 // itself, since transpose(inverse(R)) = transpose(transpose(R)) = R.

551 const double c = std::cos(0.7);

552 const double s = std::sin(0.7);

553 const fa::mat3d rot{c, -s, 0.0, s, c, 0.0, 0.0, 0.0, 1.0};

554 const fa::mat3d it = fa::inverse_transpose(rot);

555 for (std::size_t i = 0; i < 9; ++i) { EXPECT_NEAR(it.data()[i], rot.data()[i], 1e-12); }

556 // For a non-uniform scale S = diag(2, 4), inverse_transpose is diag(1/2, 1/4).

557 const fa::mat2d scale{2.0, 0.0, 0.0, 4.0};

558 const fa::mat2d nrm = fa::inverse_transpose(scale);

559 EXPECT_NEAR(nrm(0, 0), 0.5, 1e-12);

560 EXPECT_NEAR(nrm(1, 1), 0.25, 1e-12);

561 // Singular still throws (via inverse).

562 EXPECT_THROW((void)fa::inverse_transpose(fa::mat3d{}), std::domain_error);

563}

564

565// ---- Enum subscripting: a scoped enum names an axis, and only when indexing --------------------

566

567namespace {

568/// A caller's scoped enum. It stays strongly typed everywhere except at a subscript, which is the

569/// whole point of ndarray::Subscript.

570enum class Axis : std::uint8_t { X = 0, Y = 1, Z = 2 };

571enum class Basis : std::uint8_t { Right = 0, Up = 1, Forward = 2 };

572} // namespace

573

574TEST(Fixarray, EnumIndexingOnVectorsAndMatrices) {

575 fa::vec3f v{7.0F, 8.0F, 9.0F};

576 // read a component by name

577 EXPECT_FLOAT_EQ(v[Axis::X], 7.0F);

578 EXPECT_FLOAT_EQ(v[Axis::Z], 9.0F);

579 // write by name (non-const overload)

580 v[Axis::Y] = 42.0F;

581 EXPECT_FLOAT_EQ(v[1], 42.0F);

582 // const overload

583 const fa::vec3f& cv = v;

584 EXPECT_FLOAT_EQ(cv[Axis::Y], 42.0F);

585 // a plain integer still resolves the ordinary overload

586 EXPECT_FLOAT_EQ(v[std::size_t{0}], 7.0F);

587

588 fa::mat3f m = fa::mat3f::identity();

589 // both indices named

590 EXPECT_FLOAT_EQ(m(Axis::Y, Axis::Y), 1.0F);

591 EXPECT_FLOAT_EQ(m(Axis::X, Axis::Y), 0.0F);

592 // mixed: one enum, one integer

593 EXPECT_FLOAT_EQ(m(Axis::Z, std::size_t{2}), 1.0F);

594 EXPECT_FLOAT_EQ(m(std::size_t{0}, Axis::X), 1.0F);

595 // write by name

596 m(Axis::X, Axis::Z) = 5.0F;

597 EXPECT_FLOAT_EQ(m(0, 2), 5.0F);

598 // const overloads (both-enum and mixed)

599 const fa::mat3f& cm = m;

600 EXPECT_FLOAT_EQ(cm(Axis::X, Axis::Z), 5.0F);

601 EXPECT_FLOAT_EQ(cm(Axis::Y, std::size_t{1}), 1.0F);

602 EXPECT_FLOAT_EQ(cm(std::size_t{2}, Axis::Z), 1.0F);

603}

604

605TEST(Fixarray, NamedRowsAndColumns) {

606 const fa::mat3f id = fa::mat3f::identity();

607 // a basis vector by name — the axis an enum was made for

608 EXPECT_TRUE(fa::column(id, Basis::Forward) == (fa::vec3f{0.0F, 0.0F, 1.0F}));

609 EXPECT_TRUE(fa::column(id, Basis::Right) == (fa::vec3f{1.0F, 0.0F, 0.0F}));

610 // a plain integer index still works

611 EXPECT_TRUE(fa::column(id, 1) == (fa::vec3f{0.0F, 1.0F, 0.0F}));

612 EXPECT_TRUE(fa::row(id, 2) == (fa::vec3f{0.0F, 0.0F, 1.0F}));

613 EXPECT_TRUE(fa::row(id, Axis::X) == (fa::vec3f{1.0F, 0.0F, 0.0F}));

614

615 // On a real (column-major) transform, column j is the image of basis vector j.

616 const fa::mat3f t{2.0F, 0.0F, 1.0F, 0.0F, 3.0F, 2.0F, 0.0F, 0.0F, 1.0F};

617 EXPECT_TRUE(fa::column(t, Axis::X) == (fa::vec3f{2.0F, 0.0F, 0.0F})); // where x-hat lands

618 EXPECT_TRUE(fa::row(t, Axis::X) == (fa::vec3f{2.0F, 0.0F, 1.0F}));

619 // A non-square matrix: row length is the column count and vice-versa.

620 const fa::Mat<double, 2, 3> wide{1.0, 2.0, 3.0, 4.0, 5.0, 6.0};

621 EXPECT_TRUE(fa::row(wide, 1) == (fa::vec3d{4.0, 5.0, 6.0}));

622 EXPECT_TRUE(fa::column(wide, 2) == (fa::vec2d{3.0, 6.0}));

623}

624

625// ---- from_indices: the one-pass elementwise builder the component-wise ops ride on --------------

626

627TEST(Fixarray, FromIndices) {

628 // A vector built from its flat index.

629 const fa::vec4f v = fa::vec4f::from_indices([](std::size_t i) { return static_cast<float>(i * i); });

630 EXPECT_TRUE(v == (fa::vec4f{0.0F, 1.0F, 4.0F, 9.0F}));

631

632 // For a matrix the index runs over STORAGE order (column-major), so building the identity by

633 // "1 on the diagonal" means indices divisible by rows+1 — the same fact identity() uses.

634 const fa::mat3f id =

635 fa::mat3f::from_indices([](std::size_t i) { return i % 4 == 0 ? 1.0F : 0.0F; });

636 EXPECT_TRUE(id == fa::mat3f::identity());

637

638 // It is usable at compile time.

639 constexpr fa::vec3d ramp = fa::vec3d::from_indices([](std::size_t i) { return static_cast<double>(i); });

640 static_assert(ramp[2] == 2.0);

641

642 // Column-major storage is observable: element k of the buffer is what f(k) returned.

643 const fa::mat2f m = fa::mat2f::from_indices([](std::size_t i) { return static_cast<float>(i); });

644 EXPECT_EQ(m.data()[0], 0.0F);

645 EXPECT_EQ(m.data()[3], 3.0F);

646 EXPECT_EQ(m(0, 0), 0.0F);

647 EXPECT_EQ(m(0, 1), 2.0F); // flat index 2 is (row 0, col 1) in column-major

648}

649

650// ---- display: to_string and the stream operator, in the NDArray's nested-bracket form -----------

651

652TEST(Fixarray, ToStringMatchesTheNDArrayRendering) {

653 // A vector is one bracket level; elements go through the SHARED scalar formatter

654 // (1.5 prints "1.5", a whole number prints with no trailing ".0").

655 EXPECT_EQ(fa::to_string(fa::vec3f{1.5F, -2.0F, 3.0F}), "[1.5, -2, 3]");

656 // A matrix renders in reading (row, column) order regardless of the column-major storage.

657 EXPECT_EQ(fa::to_string(fa::mat2f{1.0F, 2.0F, 3.0F, 4.0F}), "[[1, 2], [3, 4]]");

658 // The double instantiation formats identically.

659 EXPECT_EQ(fa::to_string(fa::vec2d{0.25, 42.0}), "[0.25, 42]");

660}

661

662TEST(Fixarray, StreamInsertionUsesTheToStringForm) {

663 std::ostringstream vs;

664 vs << fa::vec3f{1.0F, 2.5F, -3.0F};

665 EXPECT_EQ(vs.str(), "[1, 2.5, -3]");

666 std::ostringstream ms;

667 ms << fa::mat2f{1.0F, 2.0F, 3.0F, 4.0F};

668 EXPECT_EQ(ms.str(), "[[1, 2], [3, 4]]");

669}

670

671// ---- builtins::index — what cheatah's value-position subscript v[i] / m[i, j] lowers to ---------

672

673TEST(Fixarray, BuiltinsIndexLowersSubscripts) {

674 const fa::vec3f v{7.0F, 8.0F, 9.0F};

675 EXPECT_FLOAT_EQ(cheatah::builtins::index(v, std::size_t{1}), 8.0F);

676 EXPECT_FLOAT_EQ(cheatah::builtins::index(v, Axis::Z), 9.0F); // enum labels work here too

677

678 const fa::mat2f m{1.0F, 2.0F, 3.0F, 4.0F}; // reading order

679 EXPECT_FLOAT_EQ(cheatah::builtins::index(m, std::size_t{1}, std::size_t{0}), 3.0F);

680 EXPECT_FLOAT_EQ(cheatah::builtins::index(m, Axis::X, Axis::Y), 2.0F); // (row 0, col 1)

681}