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 <sstream>

20#include <stdexcept>

21#include <type_traits>

22#include <vector>

24#include <gtest/gtest.h>

26#include "ndarray.hpp"

27#include "routines.hpp"

29namespace fa = cheatah::fixarray;

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

31namespace nd = cheatah::ndarray;

33namespace {

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

36template <class M>

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

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

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

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

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

42 return ::testing::AssertionFailure()

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

44 }

45 }

46 return ::testing::AssertionSuccess();

47}

49} // namespace

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

53TEST(Fixarray, IsAPlainValueOfExactlyItsElements) {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

79 SUCCEED();

80}

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

84TEST(Fixarray, DefaultIsZero) {

85 const fa::vec3f v;

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

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

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

89 const fa::mat2d m;

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

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

92}

94TEST(Fixarray, VectorIndexing) {

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

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

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

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

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

100 const fa::vec3f& cv = v;

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

102

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

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

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

106}

107

108TEST(Fixarray, MatrixIndexing) {

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

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

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

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

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

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

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

116 const fa::mat2f& cm = m;

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

118}

119

120TEST(Fixarray, Data) {

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

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

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

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

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

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

127

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

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

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

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

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

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

134

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

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

137 const fa::mat2f& cm = m;

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

139

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

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

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

143}

144

145TEST(Fixarray, Identity) {

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

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

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

149

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

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

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

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

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

155 }

156 }

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

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

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

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

161}

162

163TEST(Fixarray, Filled) {

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

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

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

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

168}

169

170TEST(Fixarray, Equality) {

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

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

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

174 EXPECT_TRUE(a == b);

175 EXPECT_FALSE(a == c);

176 EXPECT_TRUE(a != c);

177 EXPECT_FALSE(a != b);

178}

179

180// ---- Arithmetic ---------------------------------------------------------------------------------

181

182TEST(Fixarray, Arithmetic) {

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

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

185

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

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

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

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

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

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

192

193 fa::vec3f m = a;

194 m += b;

195 EXPECT_TRUE(m == a + b);

196 m -= b;

197 EXPECT_TRUE(m == a);

198 m *= 3.0F;

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

200 m /= 3.0F;

201 EXPECT_TRUE(m == a);

202

203 // Doubles behave the same.

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

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

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

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

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

209 dm *= 2.0;

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

211 dm /= 2.0;

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

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

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

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

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

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

218

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

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

221 small /= 2.0;

222 EXPECT_EQ(small[1], 2.0);

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

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

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

226}

227

228// ---- Vector products ----------------------------------------------------------------------------

229

230TEST(Fixarray, DotAndCross) {

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

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

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

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

235

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

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

238

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

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

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

242 // ...and anticommutative.

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

244

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

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

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

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

249 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);

250}

251

252TEST(Fixarray, NormAndNormalize) {

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

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

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

256

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

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

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

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

261

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

263

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

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

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

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

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

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

270

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

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

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

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

275}

276

277// ---- Matrix products, transpose, trace ----------------------------------------------------------

278

279TEST(Fixarray, Matmul) {

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

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

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

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

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

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

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

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

288

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

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

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

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

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

294

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

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

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

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

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

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

301

302 // Matrix times vector.

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

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

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

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

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

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

309 EXPECT_DOUBLE_EQ(rect[0], 6.0);

310 EXPECT_DOUBLE_EQ(rect[1], 15.0);

311}

312

313TEST(Fixarray, TransposeAndTrace) {

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

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

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

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

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

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

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

321

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

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

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

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

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

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

328}

329

330// ---- Determinant and inverse --------------------------------------------------------------------

331

332TEST(Fixarray, DeterminantAndInverse) {

333 // 2x2

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

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

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

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

338

339 // 3x3

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

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

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

343

344 // 4x4

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

346 2.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0};

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

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

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

350

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

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

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

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

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

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

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

358 2.0F, 0.0F, 1.0F, 1.0F, 1.0F, 1.0F, 1.0F, 2.0F}),

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

360

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

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

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

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

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

366

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

368 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}),

369 std::domain_error);

370}

371

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

373

374TEST(Fixarray, WorksBeyondTheAliasedSizes) {

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

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

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

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

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

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

381 b[i] = 1.0;

382 }

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

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

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

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

387

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

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

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

391

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

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

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

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

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

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

398 EXPECT_TRUE(through == five);

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

400}

401

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

403

404TEST(Fixarray, AgreesWithTheDynamicNDArray) {

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

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

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

408

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

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

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

412

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

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

415

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

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

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

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

420 }

421 }

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

423}

424

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

426

427TEST(Fixarray, Geometry) {

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

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

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

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

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

433

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

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

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

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

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

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

440

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

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

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

444 // Total internal reflection returns the zero vector.

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

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

447

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

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

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

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

452}

453

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

455

456TEST(Fixarray, CommonUnary) {

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

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

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

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

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

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

463}

464

465TEST(Fixarray, MinMaxClamp) {

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

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

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

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

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

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

472

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

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

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

476 // per-element bounds

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

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

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

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

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

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

483}

484

485TEST(Fixarray, MixStep) {

486 // mix with a scalar factor is a lerp

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

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

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

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

491 // per-element factor

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

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

494

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

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

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

498

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

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

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

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

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

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

505 // monotone and within [0,1]

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

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

508 EXPECT_GE(q[0], 0.0);

509 EXPECT_LE(q[1], 1.0);

510}

511

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

513

514TEST(Fixarray, MatrixExtras) {

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

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

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

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

519

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

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

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

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

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

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

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

527

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

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

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

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

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

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

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

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

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

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

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

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

540 // Singular still throws (via inverse).

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

542}

543

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

545

546namespace {

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

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

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

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

551} // namespace

552

553TEST(Fixarray, EnumIndexingOnVectorsAndMatrices) {

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

555 // read a component by name

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

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

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

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

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

561 // const overload

562 const fa::vec3f& cv = v;

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

564 // a plain integer still resolves the ordinary overload

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

566

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

568 // both indices named

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

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

571 // mixed: one enum, one integer

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

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

574 // write by name

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

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

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

578 const fa::mat3f& cm = m;

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

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

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

582}

583

584TEST(Fixarray, NamedRowsAndColumns) {

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

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

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

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

589 // a plain integer index still works

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

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

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

593

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

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

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

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

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

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

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

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

602}

603

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

605

606TEST(Fixarray, FromIndices) {

607 // A vector built from its flat index.

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

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

610

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

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

613 const fa::mat3f id =

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

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

616

617 // It is usable at compile time.

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

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

620

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

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

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

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

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

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

627}

628

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

630

631TEST(Fixarray, ToStringMatchesTheNDArrayRendering) {

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

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

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

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

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

637 // The double instantiation formats identically.

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

639}

640

641TEST(Fixarray, StreamInsertionUsesTheToStringForm) {

642 std::ostringstream vs;

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

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

645 std::ostringstream ms;

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

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

648}

649

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

651

652TEST(Fixarray, BuiltinsIndexLowersSubscripts) {

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

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

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

656

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

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

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

660}