1/* Template (no include guard: included once per instantiation; see pathtrace/real.h): programs/
2 * pathtrace's vector math, sqrt and cis over the instantiation's float operations. */
3#include "pathtrace/real.h"
4
5#ifndef PT_API
6#define PT_API static inline
7#endif
8
9static inline pt_v PT_FN(v3)(float x, float y, float z) {
10    pt_v v;
11    v.x = x;
12    v.y = y;
13    v.z = z;
14    return v;
15}
16static inline float PT_FN(ldr)(const uint32_t *sc, size_t i) { return PT_K(sc[i]); }
17static inline pt_v PT_FN(ldv)(const uint32_t *sc, size_t i) {
18    return PT_FN(v3)(PT_K(sc[i]), PT_K(sc[i + 1]), PT_K(sc[i + 2]));
19}
20static inline pt_v PT_FN(vadd)(pt_v a, pt_v b) { return PT_FN(v3)(PT_ADD(a.x, b.x), PT_ADD(a.y, b.y), PT_ADD(a.z, b.z)); }
21static inline pt_v PT_FN(vsub)(pt_v a, pt_v b) { return PT_FN(v3)(PT_SUB(a.x, b.x), PT_SUB(a.y, b.y), PT_SUB(a.z, b.z)); }
22static inline pt_v PT_FN(vmul)(pt_v a, pt_v b) { return PT_FN(v3)(PT_MUL(a.x, b.x), PT_MUL(a.y, b.y), PT_MUL(a.z, b.z)); }
23static inline pt_v PT_FN(vscale)(pt_v a, float s) { return PT_FN(v3)(PT_MUL(a.x, s), PT_MUL(a.y, s), PT_MUL(a.z, s)); }
24static inline pt_v PT_FN(vneg)(pt_v a) { return PT_FN(v3)(PT_NEG(a.x), PT_NEG(a.y), PT_NEG(a.z)); }
25/* (a.x b.x + a.y b.y) + a.z b.z: never a builtin dot (which may fuse or reorder). */
26static inline float PT_FN(dot)(pt_v a, pt_v b) {
27    return PT_ADD(PT_ADD(PT_MUL(a.x, b.x), PT_MUL(a.y, b.y)), PT_MUL(a.z, b.z));
28}
29static inline pt_v PT_FN(cross)(pt_v a, pt_v b) {
30    return PT_FN(v3)(PT_SUB(PT_MUL(a.y, b.z), PT_MUL(a.z, b.y)), PT_SUB(PT_MUL(a.z, b.x), PT_MUL(a.x, b.z)),
31                     PT_SUB(PT_MUL(a.x, b.y), PT_MUL(a.y, b.x)));
32}
33static inline pt_v PT_FN(zero3)(void) { return PT_FN(v3)(PT_K(PT_C_ZERO), PT_K(PT_C_ZERO), PT_K(PT_C_ZERO)); }
34
35/* sqrt by Newton from an integer first guess: three steps of y = (y + x / y) / 2 (about 1 ulp; the
36 * same bits on every backend). x <= 0 or NaN gives 0. */
37PT_API float PT_FN(sqrt)(float x) {
38    if (!PT_GT(x, PT_K(PT_C_ZERO))) return PT_K(PT_C_ZERO);
39    float y = PT_K(0x1fbd1df5u + (pt_to_bits(x) >> 1));
40    float h = PT_K(PT_C_HALF);
41    y = PT_MUL(h, PT_ADD(y, PT_DIV(x, y)));
42    y = PT_MUL(h, PT_ADD(y, PT_DIV(x, y)));
43    y = PT_MUL(h, PT_ADD(y, PT_DIV(x, y)));
44    return y;
45}
46
47static inline pt_v PT_FN(normalize)(pt_v a) {
48    float inv = PT_DIV(PT_K(PT_C_ONE), PT_FN(sqrt)(PT_FN(dot)(a, a)));
49    return PT_FN(vscale)(a, inv);
50}
51
52/* (cos, sin) of a quarter turn's quadrant h >> 30 plus the angle of h's next 24 bits within it:
53 * 2 pi h / 2^32, to 2^-26 of a turn. */
54PT_API void PT_FN(cis)(uint32_t h, float *cos_out, float *sin_out) {
55    float t = PT_MUL(PT_FROM_U32((h >> 6) & 0x00ffffffu), PT_K(PT_C_U24));
56    float a = PT_MUL(t, PT_K(PT_C_HALF_PI));
57    float a2 = PT_MUL(a, a);
58    /* The Rust's nested Horner expressions, innermost first (the same operations in the same order). */
59    float p = PT_ADD(PT_K(pt_c_sin[3]), PT_MUL(a2, PT_K(pt_c_sin[4])));
60    p = PT_ADD(PT_K(pt_c_sin[2]), PT_MUL(a2, p));
61    p = PT_ADD(PT_K(pt_c_sin[1]), PT_MUL(a2, p));
62    p = PT_ADD(PT_K(pt_c_sin[0]), PT_MUL(a2, p));
63    float s = PT_MUL(a, PT_ADD(PT_K(PT_C_ONE), PT_MUL(a2, p)));
64    float q = PT_ADD(PT_K(pt_c_cos[4]), PT_MUL(a2, PT_K(pt_c_cos[5])));
65    q = PT_ADD(PT_K(pt_c_cos[3]), PT_MUL(a2, q));
66    q = PT_ADD(PT_K(pt_c_cos[2]), PT_MUL(a2, q));
67    q = PT_ADD(PT_K(pt_c_cos[1]), PT_MUL(a2, q));
68    q = PT_ADD(PT_K(pt_c_cos[0]), PT_MUL(a2, q));
69    float c = PT_ADD(PT_K(PT_C_ONE), PT_MUL(a2, q));
70    switch (h >> 30) {
71    case 0: *cos_out = c, *sin_out = s; break;
72    case 1: *cos_out = PT_NEG(s), *sin_out = c; break;
73    case 2: *cos_out = PT_NEG(c), *sin_out = PT_NEG(s); break;
74    default: *cos_out = s, *sin_out = PT_NEG(c); break;
75    }
76}