Index: benchmark/vector/glm_vec2.cc
===================================================================
--- benchmark/vector/glm_vec2.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ benchmark/vector/glm_vec2.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,56 @@
+#include <glm/glm.hpp>
+#include <iostream>
+#define BENCH_N 1000000000
+#include "../bench.h"
+
+using namespace glm;
+
+
+int main(int argc, char* argv[]) {
+    tvec2<volatile float> v1(1,2);
+    tvec2<volatile float> v2(20.5f, 30.2f);
+
+    {
+        BENCH(
+            for (int i = 0; i < n; i ++) {
+                volatile float res = dot(v1, v2);
+                asm volatile ("":::"memory");
+            },
+            dot_b
+        );
+        std::cout << "dot: " << dot_b << std::endl;
+    }
+
+    {
+        BENCH(
+            for (int i = 0; i < n; i ++) {
+                volatile vec2 res = v1 + v2;
+                asm volatile ("":::"memory");
+            },
+            add_b
+        );
+        std::cout << "vector addition: " << add_b << std::endl;
+    }
+
+    {
+        BENCH(
+            for (int i = 0; i < n; i ++) {
+                volatile vec2 res = normalize((tvec2<float>)v1);
+                asm volatile ("":::"memory");
+            },
+            normalize_b
+        );
+        std::cout << "normalize: " << normalize_b << std::endl;
+    }
+
+    {
+        BENCH(
+            for (int i = 0; i < n; i ++) {
+                volatile vec2 res = reflect((vec2)v1, (vec2)v2);
+                asm volatile ("":::"memory");
+            },
+            reflect_b
+        );
+        std::cout << "reflect: " << reflect_b << std::endl;
+    }
+}
Index: benchmark/vector/vec2_notemplates.cfa
===================================================================
--- benchmark/vector/vec2_notemplates.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ benchmark/vector/vec2_notemplates.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,53 @@
+#include "../../libcfa/src/vec/vec2_f.hfa"
+/* #include <fstream.hfa> */
+#define BENCH_N 1000000000
+#include "../bench.h"
+
+int main(int argc, char* argv[]) {
+    volatile vec2 v1 = {1,2};
+    volatile vec2 v2 = {20.5, 30.2};
+
+    {
+        BENCH(
+            for (n) {
+                volatile float res = dot(v1, v2);
+                asm volatile ("":::"memory");
+            },
+            dot_b
+        );
+        printf("dot: %g\n", dot_b);
+    }
+
+    {
+        BENCH(
+            for (n) {
+                volatile vec2 res = v1 + v2;
+                asm volatile ("":::"memory");
+            },
+            add_b
+        );
+        printf("vector addition: %g\n", add_b);
+    }
+
+    {
+        BENCH(
+            for (n) {
+                volatile vec2 res = normalize(v1);
+                asm volatile ("":::"memory");
+            },
+            normalize_b
+        );
+        printf("normalize: %g\n", normalize_b);
+    }
+
+    {
+        BENCH(
+            for (n) {
+                volatile vec2 res = reflect(v1, v2);
+                asm volatile ("":::"memory");
+            },
+            reflect_b
+        );
+        printf("reflect : %g\n", reflect_b);
+    }
+}
Index: benchmark/vector/vec2_templates.cfa
===================================================================
--- benchmark/vector/vec2_templates.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ benchmark/vector/vec2_templates.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,53 @@
+#include "../../libcfa/src/vec/vec2.hfa"
+/* #include <fstream.hfa> */
+/* #define BENCH_N 1000000 */
+#include "../bench.h"
+
+int main(int argc, char* argv[]) {
+    volatile vec2(float) v1 = {1,2};
+    volatile vec2(float) v2 = {20.5f, 30.2f};
+
+    {
+        BENCH(
+            for (n) {
+                volatile float res = dot(v1, v2);
+                asm volatile ("":::"memory");
+            },
+            dot_b
+        );
+        printf("dot: %g\n", dot_b);
+    }
+
+    {
+        BENCH(
+            for (n) {
+                volatile vec2(float) res = v1 + v2;
+                asm volatile ("":::"memory");
+            },
+            add_b
+        );
+        printf("vector addition: %g\n", add_b);
+    }
+
+    {
+        BENCH(
+            for (n) {
+                volatile vec2(float) res = normalize(v1);
+                asm volatile ("":::"memory");
+            },
+            normalize_b
+        );
+        printf("normalize: %g\n", normalize_b);
+    }
+
+    {
+        BENCH(
+            for (n) {
+                volatile vec2(float) res = reflect(v1, v2);
+                asm volatile ("":::"memory");
+            },
+            reflect_b
+        );
+        printf("reflect : %g\n", reflect_b);
+    }
+}
Index: libcfa/src/Makefile.am
===================================================================
--- libcfa/src/Makefile.am	(revision 44e37ef2cacfe5adf2c7846da0808b66f7a21a4e)
+++ libcfa/src/Makefile.am	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -58,5 +58,9 @@
 	concurrency/iofwd.hfa \
 	containers/list.hfa \
-	containers/stackLockFree.hfa
+	containers/stackLockFree.hfa \
+	vec/vec.hfa \
+	vec/vec2.hfa \
+	vec/vec3.hfa \
+	vec/vec4.hfa
 
 inst_headers_src = \
Index: libcfa/src/vec/vec.hfa
===================================================================
--- libcfa/src/vec/vec.hfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ libcfa/src/vec/vec.hfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,135 @@
+//
+// Cforall Version 1.0.0 Copyright (C) 2021 University of Waterloo
+//
+// The contents of this file are covered under the licence agreement in the
+// file "LICENCE" distributed with Cforall.
+//
+// io/types.hfa --
+//
+// Author           : Dimitry Kobets
+// Created On       :
+// Last Modified By :
+// Last Modified On :
+// Update Count     :
+//
+
+#pragma once
+
+#include <math.hfa>
+
+trait fromint(otype T) {
+    void ?{}(T&, int);
+};
+trait zeroinit(otype T) {
+    void ?{}(T&, zero_t);
+};
+trait zero_assign(otype T) {
+    T ?=?(T&, zero_t);
+};
+trait subtract(otype T) {
+    T ?-?(T, T);
+};
+trait negate(otype T) {
+    T -?(T);
+};
+trait add(otype T) {
+    T ?+?(T, T);
+};
+trait multiply(otype T) {
+    T ?*?(T, T);
+};
+trait divide(otype T) {
+    T ?/?(T, T);
+};
+trait lessthan(otype T) {
+    int ?<?(T, T);
+};
+trait equality(otype T) {
+    int ?==?(T, T);
+};
+trait sqrt(otype T) {
+    T sqrt(T);
+};
+
+static inline {
+// int
+int ?=?(int& n, zero_t) { return n = 0.f; }
+// unsigned int
+int ?=?(unsigned int& n, zero_t) { return n = 0.f; }
+/* float */
+void ?{}(float& a, int b) { a = b; }
+float ?=?(float& n, zero_t) { return n = 0.f; }
+/* double */
+void ?{}(double& a, int b) { a = b; }
+double ?=?(double& n, zero_t) { return n = 0L; }
+// long double
+void ?{}(long double& a, int b) { a = b; }
+long double ?=?(long double& n, zero_t) { return n = 0L; }
+}
+
+trait dottable(otype V, otype T) {
+    T dot(V, V);
+};
+
+static inline {
+
+forall(otype T | sqrt(T), otype V | dottable(V, T))
+T length(V v) {
+   return sqrt(dot(v, v));
+}
+
+forall(otype T, otype V | dottable(V, T))
+T length_squared(V v) {
+   return dot(v, v);
+}
+
+forall(otype T, otype V | { T length(V); } | subtract(V))
+T distance(V v1, V v2) {
+    return length(v1 - v2);
+}
+
+forall(otype T, otype V | { T length(V); V ?/?(V, T); })
+V normalize(V v) {
+    return v / length(v);
+}
+
+// Project vector u onto vector v
+forall(otype T, otype V | dottable(V, T) | { V normalize(V); V ?*?(V, T); })
+V project(V u, V v) {
+    V v_norm = normalize(v);
+    return v_norm * dot(u, v_norm);
+}
+
+// Reflect incident vector v with respect to surface with normal n
+forall(otype T | fromint(T), otype V | { V project(V, V); V ?*?(T, V); V ?-?(V,V); })
+V reflect(V v, V n) {
+    return v - (T){2} * project(v, n);
+}
+
+// Refract incident vector v with respect to surface with normal n
+// eta is the ratio of indices of refraction between starting material and
+// entering material (i.e., from air to water, eta = 1/1.33)
+// v and n must already be normalized
+forall(otype T | fromint(T) | subtract(T) | multiply(T) | add(T) | lessthan(T) | sqrt(T),
+       otype V | dottable(V, T) | { V ?*?(T, V); V ?-?(V,V); void ?{}(V&, zero_t); })
+V refract(V v, V n, T eta) {
+    T dotValue = dot(n, v);
+    T k = (T){1} - eta * eta * ((T){1} - dotValue * dotValue);
+    if (k < (T){0}) {
+        return 0;
+    }
+    return eta * v - (eta * dotValue + sqrt(k)) * n;
+}
+
+// Given a perturbed normal and a geometric normal,
+// flip the perturbed normal if the geometric normal is pointing away
+// from the observer.
+// n is the perturbed vector that we want to align
+// i is the incident vector
+// ng is the geometric normal of the surface
+forall(otype T | lessthan(T) | zeroinit(T), otype V | dottable(V, T) | negate(V))
+V faceforward(V n, V i, V ng) {
+    return dot(ng, i) < (T){0} ? n : -n;
+}
+
+} // inline
Index: libcfa/src/vec/vec2.hfa
===================================================================
--- libcfa/src/vec/vec2.hfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ libcfa/src/vec/vec2.hfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,288 @@
+//
+// Cforall Version 1.0.0 Copyright (C) 2021 University of Waterloo
+//
+// The contents of this file are covered under the licence agreement in the
+// file "LICENCE" distributed with Cforall.
+//
+// io/types.hfa --
+//
+// Author           : Dimitry Kobets
+// Created On       :
+// Last Modified By :
+// Last Modified On :
+// Update Count     :
+//
+
+#pragma once
+
+#include <iostream.hfa>
+#include "vec.hfa"
+
+forall (otype T) {
+    struct vec2 {
+        T x, y;
+    };
+}
+
+forall (otype T) {
+    static inline {
+
+    void ?{}(vec2(T)& v, T x, T y) {
+        v.[x, y] = [x, y];
+    }
+
+    forall(| zero_assign(T))
+    void ?{}(vec2(T)& vec, zero_t) with (vec) {
+        x = y = 0;
+    }
+
+    void ?{}(vec2(T)& vec, T val) with (vec) {
+        x = y = val;
+    }
+
+    void ?{}(vec2(T)& vec, vec2(T) other) with (vec) {
+        [x,y] = other.[x,y];
+    }
+
+    void ?=?(vec2(T)& vec, vec2(T) other) with (vec) {
+        [x,y] = other.[x,y];
+    }
+    forall(| zero_assign(T))
+    void ?=?(vec2(T)& vec, zero_t) with (vec) {
+        x = y = 0;
+    }
+
+    // Primitive mathematical operations
+
+    // -
+    forall(| subtract(T)) {
+    vec2(T) ?-?(vec2(T) u, vec2(T) v) {
+        return [u.x - v.x, u.y - v.y];
+    }
+    vec2(T)& ?-=?(vec2(T)& u, vec2(T) v) {
+        u = u - v;
+        return u;
+    }
+    }
+    forall(| negate(T))
+    vec2(T) -?(vec2(T) v) with (v) {
+        return [-x, -y];
+    }
+
+    forall(| { T --?(T&); }) {
+    vec2(T)& --?(vec2(T)& v) {
+        --v.x;
+        --v.y;
+        return v;
+    }
+    vec2(T) ?--(vec2(T)& v) {
+        vec2(T) copy = v;
+        --v;
+        return copy;
+    }
+    }
+
+    // +
+    forall(| add(T)) {
+    vec2(T) ?+?(vec2(T) u, vec2(T) v) {
+        return [u.x + v.x, u.y + v.y];
+    }
+    vec2(T)& ?+=?(vec2(T)& u, vec2(T) v) {
+        u = u + v;
+        return u;
+    }
+    }
+
+    forall(| { T ++?(T&); }) {
+    vec2(T)& ++?(vec2(T)& v) {
+        ++v.x;
+        ++v.y;
+        return v;
+    }
+    vec2(T) ?++(vec2(T)& v) {
+        vec2(T) copy = v;
+        ++v;
+        return copy;
+    }
+    }
+
+    // *
+    forall(| multiply(T)) {
+    vec2(T) ?*?(vec2(T) v, T scalar) with (v) {
+        return [x * scalar, y * scalar];
+    }
+    vec2(T) ?*?(T scalar, vec2(T) v) {
+        return v * scalar;
+    }
+    vec2(T) ?*?(vec2(T) u, vec2(T) v) {
+        return [u.x * v.x, u.y * v.y];
+    }
+    vec2(T)& ?*=?(vec2(T)& v, T scalar) {
+        v = v * scalar;
+        return v;
+    }
+    vec2(T) ?*=?(vec2(T)& u, vec2(T) v) {
+        u = u * v;
+        return u;
+    }
+    }
+
+    // /
+    forall(| divide(T)) {
+    vec2(T) ?/?(vec2(T) v, T scalar) with (v) {
+        return [x / scalar, y / scalar];
+    }
+    vec2(T) ?/?(vec2(T) u, vec2(T) v) {
+        return [u.x / v.x, u.y / v.y];
+    }
+    vec2(T)& ?/=?(vec2(T)& v, T scalar) {
+        v = v / scalar;
+        return v;
+    }
+    vec2(T) ?/=?(vec2(T)& u, vec2(T) v) {
+        u = u / v;
+        return u;
+    }
+    }
+
+    // %
+    forall(| { T ?%?(T,T); }) {
+    vec2(T) ?%?(vec2(T) v, T scalar) with (v) {
+        return [x % scalar, y % scalar];
+    }
+    vec2(T)& ?%=?(vec2(T)& u, T scalar) {
+        u = u % scalar;
+        return u;
+    }
+    vec2(T) ?%?(vec2(T) u, vec2(T) v) {
+        return [u.x % v.x, u.y % v.y];
+    }
+    vec2(T)& ?%=?(vec2(T)& u, vec2(T) v) {
+        u = u % v;
+        return u;
+    }
+    }
+
+    // &
+    forall(| { T ?&?(T,T); }) {
+    vec2(T) ?&?(vec2(T) v, T scalar) with (v) {
+        return [x & scalar, y & scalar];
+    }
+    vec2(T)& ?&=?(vec2(T)& u, T scalar) {
+        u = u & scalar;
+        return u;
+    }
+    vec2(T) ?&?(vec2(T) u, vec2(T) v) {
+        return [u.x & v.x, u.y & v.y];
+    }
+    vec2(T)& ?&=?(vec2(T)& u, vec2(T) v) {
+        u = u & v;
+        return u;
+    }
+    }
+
+    // |
+    forall(| { T ?|?(T,T); }) {
+    vec2(T) ?|?(vec2(T) v, T scalar) with (v) {
+        return [x | scalar, y | scalar];
+    }
+    vec2(T)& ?|=?(vec2(T)& u, T scalar) {
+        u = u | scalar;
+        return u;
+    }
+    vec2(T) ?|?(vec2(T) u, vec2(T) v) {
+        return [u.x | v.x, u.y | v.y];
+    }
+    vec2(T)& ?|=?(vec2(T)& u, vec2(T) v) {
+        u = u | v;
+        return u;
+    }
+    }
+
+    // ^
+    forall(| { T ?^?(T,T); }) {
+    vec2(T) ?^?(vec2(T) v, T scalar) with (v) {
+        return [x ^ scalar, y ^ scalar];
+    }
+    vec2(T)& ?^=?(vec2(T)& u, T scalar) {
+        u = u ^ scalar;
+        return u;
+    }
+    vec2(T) ?^?(vec2(T) u, vec2(T) v) {
+        return [u.x ^ v.x, u.y ^ v.y];
+    }
+    vec2(T)& ?^=?(vec2(T)& u, vec2(T) v) {
+        u = u ^ v;
+        return u;
+    }
+    }
+
+    // <<
+    forall(| { T ?<<?(T,T); }) {
+    vec2(T) ?<<?(vec2(T) v, T scalar) with (v) {
+        return [x << scalar, y << scalar];
+    }
+    vec2(T)& ?<<=?(vec2(T)& u, T scalar) {
+        u = u << scalar;
+        return u;
+    }
+    vec2(T) ?<<?(vec2(T) u, vec2(T) v) {
+        return [u.x << v.x, u.y << v.y];
+    }
+    vec2(T)& ?<<=?(vec2(T)& u, vec2(T) v) {
+        u = u << v;
+        return u;
+    }
+    }
+
+    // >>
+    forall(| { T ?>>?(T,T); }) {
+    vec2(T) ?>>?(vec2(T) v, T scalar) with (v) {
+        return [x >> scalar, y >> scalar];
+    }
+    vec2(T)& ?>>=?(vec2(T)& u, T scalar) {
+        u = u >> scalar;
+        return u;
+    }
+    vec2(T) ?>>?(vec2(T) u, vec2(T) v) {
+        return [u.x >> v.x, u.y >> v.y];
+    }
+    vec2(T)& ?>>=?(vec2(T)& u, vec2(T) v) {
+        u = u >> v;
+        return u;
+    }
+    }
+
+    // ~
+    forall(| { T ~?(T); })
+    vec2(T) ~?(vec2(T) v) with (v) {
+        return [~v.x, ~v.y];
+    }
+
+    // relational
+    forall(| equality(T)) {
+    bool ?==?(vec2(T) u, vec2(T) v) with (u) {
+        return x == v.x && y == v.y;
+    }
+    bool ?!=?(vec2(T) u, vec2(T) v) {
+        return !(u == v);
+    }
+    }
+
+    // Geometric functions
+    forall(| add(T) | multiply(T))
+    T dot(vec2(T) u, vec2(T) v) {
+        return u.x * v.x + u.y * v.y;
+    }
+
+    } // static inline
+}
+
+forall(dtype ostype, otype T | writeable(T, ostype)) {
+    ostype & ?|?(ostype & os, vec2(T) v) with (v) {
+        return os | '<' | x | ',' | y | '>';
+    }
+    void ?|?(ostype & os, vec2(T) v ) with (v) {
+        (ostype &)(os | v); ends(os);
+    }
+}
Index: libcfa/src/vec/vec3.hfa
===================================================================
--- libcfa/src/vec/vec3.hfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ libcfa/src/vec/vec3.hfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,297 @@
+//
+// Cforall Version 1.0.0 Copyright (C) 2021 University of Waterloo
+//
+// The contents of this file are covered under the licence agreement in the
+// file "LICENCE" distributed with Cforall.
+//
+// io/types.hfa --
+//
+// Author           : Dimitry Kobets
+// Created On       :
+// Last Modified By :
+// Last Modified On :
+// Update Count     :
+//
+
+#pragma once
+
+#include <iostream.hfa>
+#include "vec.hfa"
+
+forall (otype T) {
+    struct vec3 {
+        T x, y, z;
+    };
+}
+
+forall (otype T) {
+    static inline {
+
+    void ?{}(vec3(T)& v, T x, T y, T z) {
+        v.[x, y, z] = [x, y, z];
+    }
+
+    forall(| zero_assign(T))
+    void ?{}(vec3(T)& vec, zero_t) with (vec) {
+        x = y = z = 0;
+    }
+
+    void ?{}(vec3(T)& vec, T val) with (vec) {
+        x = y = z = val;
+    }
+
+    void ?{}(vec3(T)& vec, vec3(T) other) with (vec) {
+        [x,y,z] = other.[x,y,z];
+    }
+
+    void ?=?(vec3(T)& vec, vec3(T) other) with (vec) {
+        [x,y,z] = other.[x,y,z];
+    }
+    forall(| zero_assign(T))
+    void ?=?(vec3(T)& vec, zero_t) with (vec) {
+        x = y = z = 0;
+    }
+
+    // Primitive mathematical operations
+
+    // -
+    forall(| subtract(T)) {
+    vec3(T) ?-?(vec3(T) u, vec3(T) v) {
+        return [u.x - v.x, u.y - v.y, u.z - v.z];
+    }
+    vec3(T)& ?-=?(vec3(T)& u, vec3(T) v) {
+        u = u - v;
+        return u;
+    }
+    }
+    forall(| negate(T)) {
+    vec3(T) -?(vec3(T) v) with (v) {
+        return [-x, -y, -z];
+    }
+    }
+    forall(| { T --?(T&); }) {
+    vec3(T)& --?(vec3(T)& v) {
+        --v.x;
+        --v.y;
+        --v.z;
+        return v;
+    }
+    vec3(T) ?--(vec3(T)& v) {
+        vec3(T) copy = v;
+        --v;
+        return copy;
+    }
+    }
+
+    // +
+    forall(| add(T)) {
+    vec3(T) ?+?(vec3(T) u, vec3(T) v) {
+        return [u.x + v.x, u.y + v.y, u.z + v.z];
+    }
+    vec3(T)& ?+=?(vec3(T)& u, vec3(T) v) {
+        u = u + v;
+        return u;
+    }
+    }
+
+    forall(| { T ++?(T&); }) {
+    vec3(T)& ++?(vec3(T)& v) {
+        ++v.x;
+        ++v.y;
+        ++v.z;
+        return v;
+    }
+    vec3(T) ?++(vec3(T)& v) {
+        vec3(T) copy = v;
+        ++v;
+        return copy;
+    }
+    }
+
+    // *
+    forall(| multiply(T)) {
+    vec3(T) ?*?(vec3(T) v, T scalar) with (v) {
+        return [x * scalar, y * scalar, z * scalar];
+    }
+    vec3(T) ?*?(T scalar, vec3(T) v) {
+        return v * scalar;
+    }
+    vec3(T) ?*?(vec3(T) u, vec3(T) v) {
+        return [u.x * v.x, u.y * v.y, u.z * v.z];
+    }
+    vec3(T)& ?*=?(vec3(T)& v, T scalar) {
+        v = v * scalar;
+        return v;
+    }
+    vec3(T)& ?*=?(vec3(T)& u, vec3(T) v) {
+        u = u * v;
+        return u;
+    }
+    }
+
+    // /
+    forall(| divide(T)) {
+    vec3(T) ?/?(vec3(T) v, T scalar) with (v) {
+        return [x / scalar, y / scalar, z / scalar];
+    }
+    vec3(T) ?/?(vec3(T) u, vec3(T) v) {
+        return [u.x / v.x, u.y / v.y, u.z / v.z];
+    }
+    vec3(T)& ?/=?(vec3(T)& v, T scalar) {
+        v = v / scalar;
+        return v;
+    }
+    vec3(T)& ?/=?(vec3(T)& u, vec3(T) v) {
+        u = u / v;
+        return u;
+    }
+    }
+
+    // %
+    forall(| { T ?%?(T,T); }) {
+    vec3(T) ?%?(vec3(T) v, T scalar) with (v) {
+        return [x % scalar, y % scalar, z % scalar];
+    }
+    vec3(T)& ?%=?(vec3(T)& u, T scalar) {
+        u = u % scalar;
+        return u;
+    }
+    vec3(T) ?%?(vec3(T) u, vec3(T) v) {
+        return [u.x % v.x, u.y % v.y, u.z % v.z];
+    }
+    vec3(T)& ?%=?(vec3(T)& u, vec3(T) v) {
+        u = u % v;
+        return u;
+    }
+    }
+
+    // &
+    forall(| { T ?&?(T,T); }) {
+    vec3(T) ?&?(vec3(T) v, T scalar) with (v) {
+        return [x & scalar, y & scalar, z & scalar];
+    }
+    vec3(T)& ?&=?(vec3(T)& u, T scalar) {
+        u = u & scalar;
+        return u;
+    }
+    vec3(T) ?&?(vec3(T) u, vec3(T) v) {
+        return [u.x & v.x, u.y & v.y, u.z & v.z];
+    }
+    vec3(T)& ?&=?(vec3(T)& u, vec3(T) v) {
+        u = u & v;
+        return u;
+    }
+    }
+
+    // |
+    forall(| { T ?|?(T,T); }) {
+    vec3(T) ?|?(vec3(T) v, T scalar) with (v) {
+        return [x | scalar, y | scalar, z | scalar];
+    }
+    vec3(T)& ?|=?(vec3(T)& u, T scalar) {
+        u = u | scalar;
+        return u;
+    }
+    vec3(T) ?|?(vec3(T) u, vec3(T) v) {
+        return [u.x | v.x, u.y | v.y, u.z | v.z];
+    }
+    vec3(T)& ?|=?(vec3(T)& u, vec3(T) v) {
+        u = u | v;
+        return u;
+    }
+    }
+
+    // ^
+    forall(| { T ?^?(T,T); }) {
+    vec3(T) ?^?(vec3(T) v, T scalar) with (v) {
+        return [x ^ scalar, y ^ scalar, z ^ scalar];
+    }
+    vec3(T)& ?^=?(vec3(T)& u, T scalar) {
+        u = u ^ scalar;
+        return u;
+    }
+    vec3(T) ?^?(vec3(T) u, vec3(T) v) {
+        return [u.x ^ v.x, u.y ^ v.y, u.z ^ v.z];
+    }
+    vec3(T)& ?^=?(vec3(T)& u, vec3(T) v) {
+        u = u ^ v;
+        return u;
+    }
+    }
+
+    // <<
+    forall(| { T ?<<?(T,T); }) {
+    vec3(T) ?<<?(vec3(T) v, T scalar) with (v) {
+        return [x << scalar, y << scalar, z << scalar];
+    }
+    vec3(T)& ?<<=?(vec3(T)& u, T scalar) {
+        u = u << scalar;
+        return u;
+    }
+    vec3(T) ?<<?(vec3(T) u, vec3(T) v) {
+        return [u.x << v.x, u.y << v.y, u.z << v.z];
+    }
+    vec3(T)& ?<<=?(vec3(T)& u, vec3(T) v) {
+        u = u << v;
+        return u;
+    }
+    }
+
+    // >>
+    forall(| { T ?>>?(T,T); }) {
+    vec3(T) ?>>?(vec3(T) v, T scalar) with (v) {
+        return [x >> scalar, y >> scalar, z >> scalar];
+    }
+    vec3(T)& ?>>=?(vec3(T)& u, T scalar) {
+        u = u >> scalar;
+        return u;
+    }
+    vec3(T) ?>>?(vec3(T) u, vec3(T) v) {
+        return [u.x >> v.x, u.y >> v.y, u.z >> v.z];
+    }
+    vec3(T)& ?>>=?(vec3(T)& u, vec3(T) v) {
+        u = u >> v;
+        return u;
+    }
+    }
+
+    // ~
+    forall(| { T ~?(T); })
+    vec3(T) ~?(vec3(T) v) with (v) {
+        return [~v.x, ~v.y, ~v.z];
+    }
+
+    // relational
+    forall(| equality(T)) {
+    bool ?==?(vec3(T) u, vec3(T) v) with (u) {
+        return x == v.x && y == v.y && z == v.z;
+    }
+    bool ?!=?(vec3(T) u, vec3(T) v) {
+        return !(u == v);
+    }
+    }
+
+    // Geometric functions
+    forall(| add(T) | multiply(T))
+    T dot(vec3(T) u, vec3(T) v) {
+        return u.x * v.x + u.y * v.y + u.z * v.z;
+    }
+
+    forall(| subtract(T) | multiply(T))
+    vec3(T) cross(vec3(T) u, vec3(T) v) {
+        return (vec3(T)){ u.y * v.z - v.y * u.z,
+                          u.z * v.x - v.z * u.x,
+                          u.x * v.y - v.x * u.y };
+    }
+
+    } // static inline
+}
+
+forall(dtype ostype, otype T | writeable(T, ostype)) {
+    ostype & ?|?(ostype & os, vec3(T) v) with (v) {
+        return os | '<' | x | ',' | y | ',' | z | '>';
+    }
+    void ?|?(ostype & os, vec3(T) v ) with (v) {
+        (ostype &)(os | v); ends(os);
+    }
+}
Index: libcfa/src/vec/vec4.hfa
===================================================================
--- libcfa/src/vec/vec4.hfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ libcfa/src/vec/vec4.hfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,293 @@
+//
+// Cforall Version 1.0.0 Copyright (C) 2021 University of Waterloo
+//
+// The contents of this file are covered under the licence agreement in the
+// file "LICENCE" distributed with Cforall.
+//
+// io/types.hfa --
+//
+// Author           : Dimitry Kobets
+// Created On       :
+// Last Modified By :
+// Last Modified On :
+// Update Count     :
+//
+
+#pragma once
+
+#include <iostream.hfa>
+#include "vec.hfa"
+
+forall (otype T) {
+    struct vec4 {
+        T x, y, z, w;
+    };
+}
+
+forall (otype T) {
+    static inline {
+
+    void ?{}(vec4(T)& v, T x, T y, T z, T w) {
+        v.[x, y, z, w] = [x, y, z, w];
+    }
+
+    forall(| zero_assign(T))
+    void ?{}(vec4(T)& vec, zero_t) with (vec) {
+        x = y = z = w = 0;
+    }
+
+    void ?{}(vec4(T)& vec, T val) with (vec) {
+        x = y = z = w = val;
+    }
+
+    void ?{}(vec4(T)& vec, vec4(T) other) with (vec) {
+        [x,y,z,w] = other.[x,y,z,w];
+    }
+
+    void ?=?(vec4(T)& vec, vec4(T) other) with (vec) {
+        [x,y,z,w] = other.[x,y,z,w];
+    }
+    forall(| zero_assign(T))
+    void ?=?(vec4(T)& vec, zero_t) with (vec) {
+        x = y = z = w = 0;
+    }
+
+    // Primitive mathematical operations
+
+    // -
+    forall(| subtract(T)) {
+    vec4(T) ?-?(vec4(T) u, vec4(T) v) {
+        return [u.x - v.x, u.y - v.y, u.z - v.z, u.w - v.w];
+    }
+    vec4(T)& ?-=?(vec4(T)& u, vec4(T) v) {
+        u = u - v;
+        return u;
+    }
+    }
+    forall(| negate(T)) {
+    vec4(T) -?(vec4(T) v) with (v) {
+        return [-x, -y, -z, -w];
+    }
+    }
+    forall(| { T --?(T&); }) {
+    vec4(T)& --?(vec4(T)& v) {
+        --v.x;
+        --v.y;
+        --v.z;
+        --v.w;
+        return v;
+    }
+    vec4(T) ?--(vec4(T)& v) {
+        vec4(T) copy = v;
+        --v;
+        return copy;
+    }
+    }
+
+    // +
+    forall(| add(T)) {
+    vec4(T) ?+?(vec4(T) u, vec4(T) v) {
+        return [u.x + v.x, u.y + v.y, u.z + v.z, u.w + v.w];
+    }
+    vec4(T)& ?+=?(vec4(T)& u, vec4(T) v) {
+        u = u + v;
+        return u;
+    }
+    }
+
+    forall(| { T ++?(T&); }) {
+    vec4(T)& ++?(vec4(T)& v) {
+        ++v.x;
+        ++v.y;
+        ++v.z;
+        ++v.w;
+        return v;
+    }
+    vec4(T) ?++(vec4(T)& v) {
+        vec4(T) copy = v;
+        ++v;
+        return copy;
+    }
+    }
+
+    // *
+    forall(| multiply(T)) {
+    vec4(T) ?*?(vec4(T) v, T scalar) with (v) {
+        return [x * scalar, y * scalar, z * scalar, w * scalar];
+    }
+    vec4(T) ?*?(T scalar, vec4(T) v) {
+        return v * scalar;
+    }
+    vec4(T) ?*?(vec4(T) u, vec4(T) v) {
+        return [u.x * v.x, u.y * v.y, u.z * v.z, u.w * v.w];
+    }
+    vec4(T)& ?*=?(vec4(T)& v, T scalar) {
+        v = v * scalar;
+        return v;
+    }
+    vec4(T)& ?*=?(vec4(T)& u, vec4(T) v) {
+        u = u * v;
+        return u;
+    }
+    }
+
+    // /
+    forall(| divide(T)) {
+    vec4(T) ?/?(vec4(T) v, T scalar) with (v) {
+        return [x / scalar, y / scalar, z / scalar, w / scalar];
+    }
+    vec4(T) ?/?(vec4(T) u, vec4(T) v) {
+        return [u.x / v.x, u.y / v.y, u.z / v.z, u.w / v.w];
+    }
+    vec4(T)& ?/=?(vec4(T)& v, T scalar) {
+        v = v / scalar;
+        return v;
+    }
+    vec4(T)& ?/=?(vec4(T)& u, vec4(T) v) {
+        u = u / v;
+        return u;
+    }
+    }
+
+    // %
+    forall(| { T ?%?(T,T); }) {
+    vec4(T) ?%?(vec4(T) v, T scalar) with (v) {
+        return [x % scalar, y % scalar, z % scalar, w % scalar];
+    }
+    vec4(T)& ?%=?(vec4(T)& u, T scalar) {
+        u = u % scalar;
+        return u;
+    }
+    vec4(T) ?%?(vec4(T) u, vec4(T) v) {
+        return [u.x % v.x, u.y % v.y, u.z % v.z, u.w % v.w];
+    }
+    vec4(T)& ?%=?(vec4(T)& u, vec4(T) v) {
+        u = u % v;
+        return u;
+    }
+    }
+
+    // &
+    forall(| { T ?&?(T,T); }) {
+    vec4(T) ?&?(vec4(T) v, T scalar) with (v) {
+        return [x & scalar, y & scalar, z & scalar, w & scalar];
+    }
+    vec4(T)& ?&=?(vec4(T)& u, T scalar) {
+        u = u & scalar;
+        return u;
+    }
+    vec4(T) ?&?(vec4(T) u, vec4(T) v) {
+        return [u.x & v.x, u.y & v.y, u.z & v.z, u.w & v.w];
+    }
+    vec4(T)& ?&=?(vec4(T)& u, vec4(T) v) {
+        u = u & v;
+        return u;
+    }
+    }
+
+    // |
+    forall(| { T ?|?(T,T); }) {
+    vec4(T) ?|?(vec4(T) v, T scalar) with (v) {
+        return [x | scalar, y | scalar, z | scalar, w | scalar];
+    }
+    vec4(T)& ?|=?(vec4(T)& u, T scalar) {
+        u = u | scalar;
+        return u;
+    }
+    vec4(T) ?|?(vec4(T) u, vec4(T) v) {
+        return [u.x | v.x, u.y | v.y, u.z | v.z, u.w | v.w];
+    }
+    vec4(T)& ?|=?(vec4(T)& u, vec4(T) v) {
+        u = u | v;
+        return u;
+    }
+    }
+
+    // ^
+    forall(| { T ?^?(T,T); }) {
+    vec4(T) ?^?(vec4(T) v, T scalar) with (v) {
+        return [x ^ scalar, y ^ scalar, z ^ scalar, w ^ scalar];
+    }
+    vec4(T)& ?^=?(vec4(T)& u, T scalar) {
+        u = u ^ scalar;
+        return u;
+    }
+    vec4(T) ?^?(vec4(T) u, vec4(T) v) {
+        return [u.x ^ v.x, u.y ^ v.y, u.z ^ v.z, u.w ^ v.w];
+    }
+    vec4(T)& ?^=?(vec4(T)& u, vec4(T) v) {
+        u = u ^ v;
+        return u;
+    }
+    }
+
+    // <<
+    forall(| { T ?<<?(T,T); }) {
+    vec4(T) ?<<?(vec4(T) v, T scalar) with (v) {
+        return [x << scalar, y << scalar, z << scalar, w << scalar];
+    }
+    vec4(T)& ?<<=?(vec4(T)& u, T scalar) {
+        u = u << scalar;
+        return u;
+    }
+    vec4(T) ?<<?(vec4(T) u, vec4(T) v) {
+        return [u.x << v.x, u.y << v.y, u.z << v.z, u.w << v.w];
+    }
+    vec4(T)& ?<<=?(vec4(T)& u, vec4(T) v) {
+        u = u << v;
+        return u;
+    }
+    }
+
+    // >>
+    forall(| { T ?>>?(T,T); }) {
+    vec4(T) ?>>?(vec4(T) v, T scalar) with (v) {
+        return [x >> scalar, y >> scalar, z >> scalar, w >> scalar];
+    }
+    vec4(T)& ?>>=?(vec4(T)& u, T scalar) {
+        u = u >> scalar;
+        return u;
+    }
+    vec4(T) ?>>?(vec4(T) u, vec4(T) v) {
+        return [u.x >> v.x, u.y >> v.y, u.z >> v.z, u.w >> v.w];
+    }
+    vec4(T)& ?>>=?(vec4(T)& u, vec4(T) v) {
+        u = u >> v;
+        return u;
+    }
+    }
+
+    // ~
+    forall(| { T ~?(T); })
+    vec4(T) ~?(vec4(T) v) with (v) {
+        return [~x, ~y, ~z, ~w];
+    }
+
+    // relational
+    forall(| equality(T)) {
+    bool ?==?(vec4(T) u, vec4(T) v) with (u) {
+        return x == v.x && y == v.y && z == v.z && w == v.w;
+    }
+    bool ?!=?(vec4(T) u, vec4(T) v) {
+        return !(u == v);
+    }
+    }
+
+    // Geometric functions
+    forall(| add(T) | multiply(T))
+    T dot(vec4(T) u, vec4(T) v) {
+        return u.x * v.x + u.y * v.y + u.z * v.z + u.w * v.w;
+    }
+
+    } // static inline
+}
+
+forall(dtype ostype, otype T | writeable(T, ostype)) {
+    ostype & ?|?(ostype & os, vec4(T) v) with (v) {
+        return os | '<' | x | ',' | y | ',' | z | ',' | w | '>';
+    }
+    void ?|?(ostype & os, vec4(T) v ) with (v) {
+        (ostype &)(os | v); ends(os);
+    }
+}
+
Index: src/ResolvExpr/ResolveAssertions.cc
===================================================================
--- src/ResolvExpr/ResolveAssertions.cc	(revision 44e37ef2cacfe5adf2c7846da0808b66f7a21a4e)
+++ src/ResolvExpr/ResolveAssertions.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -397,5 +397,5 @@
 
 	/// Limit to depth of recursion of assertion satisfaction
-	static const int recursionLimit = 4;
+	static const int recursionLimit = 7;
 	/// Maximum number of simultaneously-deferred assertions to attempt concurrent satisfaction of
 	static const int deferLimit = 10;
Index: src/ResolvExpr/SatisfyAssertions.cpp
===================================================================
--- src/ResolvExpr/SatisfyAssertions.cpp	(revision 44e37ef2cacfe5adf2c7846da0808b66f7a21a4e)
+++ src/ResolvExpr/SatisfyAssertions.cpp	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -57,8 +57,8 @@
 		ast::UniqueId resnSlot;          ///< Slot for any recursive assertion IDs
 
-		AssnCandidate( 
-			const ast::SymbolTable::IdData c, const ast::Type * at, ast::TypeEnvironment && e, 
+		AssnCandidate(
+			const ast::SymbolTable::IdData c, const ast::Type * at, ast::TypeEnvironment && e,
 			ast::AssertionSet && h, ast::AssertionSet && n, ast::OpenVarSet && o, ast::UniqueId rs )
-		: cdata( c ), adjType( at ), env( std::move( e ) ), have( std::move( h ) ), 
+		: cdata( c ), adjType( at ), env( std::move( e ) ), have( std::move( h ) ),
 		  need( std::move( n ) ), open( std::move( o ) ), resnSlot( rs ) {}
 	};
@@ -73,6 +73,6 @@
 		const AssnCandidate & match;
 	};
-	
-	/// Wrapper for the deferred items from a single assertion satisfaction. 
+
+	/// Wrapper for the deferred items from a single assertion satisfaction.
 	/// Acts like an indexed list of DeferRef
 	struct DeferItem {
@@ -81,5 +81,5 @@
 		AssnCandidateList matches;
 
-		DeferItem( 
+		DeferItem(
 			const ast::VariableExpr * d, const ast::AssertionSetValue & i, AssnCandidateList && ms )
 		: expr( d ), info( i ), matches( std::move( ms ) ) {}
@@ -117,18 +117,18 @@
 		/// Initial satisfaction state for a candidate
 		SatState( CandidateRef & c, const ast::SymbolTable & syms )
-		: cand( c ), need(), newNeed(), deferred(), inferred(), costs{ Cost::zero }, 
+		: cand( c ), need(), newNeed(), deferred(), inferred(), costs{ Cost::zero },
 		  symtab( syms ) { need.swap( c->need ); }
-		
+
 		/// Update satisfaction state for next step after previous state
 		SatState( SatState && o, IterateFlag )
-		: cand( std::move( o.cand ) ), need( o.newNeed.begin(), o.newNeed.end() ), newNeed(), 
-		  deferred(), inferred( std::move( o.inferred ) ), costs( std::move( o.costs ) ), 
+		: cand( std::move( o.cand ) ), need( o.newNeed.begin(), o.newNeed.end() ), newNeed(),
+		  deferred(), inferred( std::move( o.inferred ) ), costs( std::move( o.costs ) ),
 		  symtab( o.symtab ) { costs.emplace_back( Cost::zero ); }
-		
+
 		/// Field-wise next step constructor
 		SatState(
-			CandidateRef && c, ast::AssertionSet && nn, InferCache && i, CostVec && cs, 
+			CandidateRef && c, ast::AssertionSet && nn, InferCache && i, CostVec && cs,
 			ast::SymbolTable && syms )
-		: cand( std::move( c ) ), need( nn.begin(), nn.end() ), newNeed(), deferred(), 
+		: cand( std::move( c ) ), need( nn.begin(), nn.end() ), newNeed(), deferred(),
 		  inferred( std::move( i ) ), costs( std::move( cs ) ), symtab( std::move( syms ) )
 		  { costs.emplace_back( Cost::zero ); }
@@ -143,12 +143,12 @@
 
 	/// Binds a single assertion, updating satisfaction state
-	void bindAssertion( 
-		const ast::VariableExpr * expr, const ast::AssertionSetValue & info, CandidateRef & cand, 
+	void bindAssertion(
+		const ast::VariableExpr * expr, const ast::AssertionSetValue & info, CandidateRef & cand,
 		AssnCandidate & match, InferCache & inferred
 	) {
 		const ast::DeclWithType * candidate = match.cdata.id;
-		assertf( candidate->uniqueId, 
+		assertf( candidate->uniqueId,
 			"Assertion candidate does not have a unique ID: %s", toString( candidate ).c_str() );
-		
+
 		ast::Expr * varExpr = match.cdata.combine( cand->expr->location, cand->cvtCost );
 		varExpr->result = match.adjType;
@@ -201,5 +201,5 @@
 			ast::OpenVarSet newOpen{ sat.cand->open };
 			ast::ptr< ast::Type > toType = assn.first->result;
-			ast::ptr< ast::Type > adjType = 
+			ast::ptr< ast::Type > adjType =
 				renameTyVars( adjustExprType( candidate->get_type(), newEnv, sat.symtab ), GEN_USAGE, false );
 
@@ -213,5 +213,5 @@
 				}
 
-				matches.emplace_back( 
+				matches.emplace_back(
 					cdata, adjType, std::move( newEnv ), std::move( have ), std::move( newNeed ),
 					std::move( newOpen ), crntResnSlot );
@@ -287,9 +287,9 @@
 	};
 
-	/// Replace ResnSlots with InferParams and add alternative to output list, if it meets pruning 
+	/// Replace ResnSlots with InferParams and add alternative to output list, if it meets pruning
 	/// threshold.
-	void finalizeAssertions( 
-		CandidateRef & cand, InferCache & inferred, PruneMap & thresholds, CostVec && costs, 
-		CandidateList & out 
+	void finalizeAssertions(
+		CandidateRef & cand, InferCache & inferred, PruneMap & thresholds, CostVec && costs,
+		CandidateList & out
 	) {
 		// prune if cheaper alternative for same key has already been generated
@@ -308,6 +308,6 @@
 	}
 
-	/// Combo iterator that combines candidates into an output list, merging their environments. 
-	/// Rejects an appended candidate if environments cannot be merged. See `Common/FilterCombos.h` 
+	/// Combo iterator that combines candidates into an output list, merging their environments.
+	/// Rejects an appended candidate if environments cannot be merged. See `Common/FilterCombos.h`
 	/// for description of "combo iterator".
 	class CandidateEnvMerger {
@@ -390,6 +390,6 @@
 } // anonymous namespace
 
-void satisfyAssertions( 
-	CandidateRef & cand, const ast::SymbolTable & symtab, CandidateList & out, 
+void satisfyAssertions(
+	CandidateRef & cand, const ast::SymbolTable & symtab, CandidateList & out,
 	std::vector<std::string> & errors
 ) {
@@ -418,5 +418,5 @@
 
 			// should a limit be imposed? worst case here is O(n^2) but very unlikely to happen.
-			for (unsigned resetCount = 0; ; ++resetCount) { 
+			for (unsigned resetCount = 0; ; ++resetCount) {
 				ast::AssertionList next;
 				resetTyVarRenaming();
@@ -449,5 +449,5 @@
 				// either add successful match or push back next state
 				if ( sat.newNeed.empty() ) {
-					finalizeAssertions( 
+					finalizeAssertions(
 						sat.cand, sat.inferred, thresholds, std::move( sat.costs ), out );
 				} else {
@@ -471,5 +471,5 @@
 				std::vector< CandidateEnvMerger::OutType > compatible = filterCombos(
 					sat.deferred, CandidateEnvMerger{ sat.cand->env, sat.cand->open, sat.symtab } );
-				
+
 				// fail early if no mutually-compatible assertion satisfaction
 				if ( compatible.empty() ) {
@@ -494,5 +494,5 @@
 					// set up next satisfaction state
 					CandidateRef nextCand = std::make_shared<Candidate>(
-						sat.cand->expr, std::move( compat.env ), std::move( compat.open ), 
+						sat.cand->expr, std::move( compat.env ), std::move( compat.open ),
 						ast::AssertionSet{} /* need moved into satisfaction state */,
 						sat.cand->cost, sat.cand->cvtCost );
@@ -500,8 +500,8 @@
 					ast::AssertionSet nextNewNeed{ sat.newNeed };
 					InferCache nextInferred{ sat.inferred };
-					
+
 					CostVec nextCosts{ sat.costs };
 					nextCosts.back() += compat.cost;
-								
+
 					ast::SymbolTable nextSymtab{ sat.symtab };
 
@@ -517,10 +517,10 @@
 					// either add successful match or push back next state
 					if ( nextNewNeed.empty() ) {
-						finalizeAssertions( 
+						finalizeAssertions(
 							nextCand, nextInferred, thresholds, std::move( nextCosts ), out );
 					} else {
-						nextSats.emplace_back( 
-							std::move( nextCand ), std::move( nextNewNeed ), 
-							std::move( nextInferred ), std::move( nextCosts ), 
+						nextSats.emplace_back(
+							std::move( nextCand ), std::move( nextNewNeed ),
+							std::move( nextInferred ), std::move( nextCosts ),
 							std::move( nextSymtab ) );
 					}
@@ -534,5 +534,5 @@
 		nextSats.clear();
 	}
-	
+
 	// exceeded recursion limit if reaches here
 	if ( out.empty() ) {
Index: tests/vector_math/.expect/vec2_double.txt
===================================================================
--- tests/vector_math/.expect/vec2_double.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/.expect/vec2_double.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,2 @@
+ctor(x,y):<1.79769313486232e+308,2.2250738585072e-308>
+dot(v2,v3):324.670709736639
Index: tests/vector_math/.expect/vec2_float.txt
===================================================================
--- tests/vector_math/.expect/vec2_float.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/.expect/vec2_float.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,41 @@
+ctor(x,y):<1.,2.>
+copy ctor:<1.,2.>
+assignment:<3.,4.2>
+move assignment:<1.,2.>
+zero-init:<0.,0.>
+zero-assign:<0.,0.>
+fill-ctor:<1.23,1.23>
+?-?:<0.02,0.43>
+?-=?:<3.23,-6.57>
+-?:<-3.23,6.57>
+?+?:<2.3,2.45>
+?+=?:<2.3,2.45>
+v*s:<4.5,8.25>
+s*v:<4.5,8.25>
+?*=?:<4.5,8.25>
+?*?(vec)<2.64,-17.4>
+?*=?(vec)<2.64,-17.4>
+?/?:<0.666667,-0.0333333>
+?/=?:<0.666667,-0.0333333>
+?/?(vec)<0.545455,-1.93333>
+?/=?(vec)<0.545455,-1.93333>
+++?<2.2,-2.5>
+--?<2.2,-2.5>
+?++(part 1)<1.2,-3.5>
+?++(part 2)<2.2,-2.5>
+?--(part 1)<1.2,-3.5>
+?--(part 2)<0.2,-4.5>
+dot_1:0.
+dot_2:8.6
+length:5.
+length_squared:25.
+distance:6.5146
+normalize:<0.882353,-0.470588>
+normalize_2:<1.,0.>
+project:<5.,0.>
+project_2:<5.5,5.5>
+reflect:<-5.,6.>
+refract:<5.,6.>
+refract:<3.7594,4.5592>
+faceforward_nochange:<4.,5.5>
+faceforward_flip:<-4.,-5.5>
Index: tests/vector_math/.expect/vec2_int.txt
===================================================================
--- tests/vector_math/.expect/vec2_int.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/.expect/vec2_int.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,57 @@
+ctor(x,y):<1,2>
+copy ctor:<1,2>
+assignment:<3,4>
+move assignment:<1,2>
+zero-init:<0,0>
+zero-assign:<0,0>
+fill-ctor:<123,123>
+?-?:<-1002,17>
+?-=?:<-1002,17>
+-?:<1002,-17>
+?+?:<1004,-11>
+?+=?:<1004,-11>
+v*s:<45,825>
+s*v:<45,825>
+?*=?:<45,825>
+?*?(vec)<2,-15>
+?*=?(vec)<2,-15>
+?/?:<7,-3>
+?/=?:<7,-3>
+?/?(vec)<5,-1>
+?/=?(vec)<5,-1>
+++?<2,-2>
+--?<2,-2>
+?++(part 1)<1,-3>
+?++(part 2)<2,-2>
+?--(part 1)<1,-3>
+?--(part 2)<0,-4>
+?%?(int)<0,1>
+?%=?(int)<0,1>
+?%?(vec)<-2,3>
+?%=?(vec)<-2,3>
+?&?(int)<0,8>
+?&=?(int)<0,8>
+?&?(vec)<1,6>
+?&=?(vec)<1,6>
+?|?(int)<7,3>
+?|=?(int)<7,3>
+?|?(vec)<7,1>
+?|=?(vec)<7,1>
+?^?(int)<12,0>
+?^=?(int)<12,0>
+?^?(vec)<12,5>
+?^=?(vec)<12,5>
+?<<?(int)<4,2>
+?<<=?(int)<4,2>
+?<<?(vec)<32,4>
+?<<=?(vec)<32,4>
+?>>?(int)<1,0>
+?>>=?(int)<1,0>
+?>>?(vec)<2,1>
+?>>=?(vec)<2,1>
+~?<-6,-3>
+dot_1:0
+dot_2:32
+length_squared:25
+faceforward_nochange:<4,5>
+faceforward_flip:<-4,-5>
Index: tests/vector_math/.expect/vec2_ldouble.txt
===================================================================
--- tests/vector_math/.expect/vec2_ldouble.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/.expect/vec2_ldouble.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,2 @@
+ctor(x,y):<1.18973149535723177e+4932,3.36210314311209351e-4932>
+dot(v2,v3):324.670709736639038
Index: tests/vector_math/.expect/vec2_uint.txt
===================================================================
--- tests/vector_math/.expect/vec2_uint.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/.expect/vec2_uint.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,2 @@
+ctor(x,y):<4294967295,0>
+dot(v2,v3):352
Index: tests/vector_math/.expect/vec3_float.txt
===================================================================
--- tests/vector_math/.expect/vec3_float.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/.expect/vec3_float.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,35 @@
+ctor(x,y):<1.,2.,3.>
+copy ctor:<1.,2.,3.>
+assignment:<3.,4.2,-2.>
+move assignment:<1.,2.,3.>
+zero-init:<0.,0.,0.>
+zero-assign:<0.,0.,0.>
+fill-ctor:<1.23,1.23,1.23>
+?-?:<0.02,0.43,-0.999998>
+?-=?:<0.,0.,0.>
+-?:<-0.,-0.,-0.>
+?+?:<2.3,2.45,-9.2>
+?+=?:<2.3,2.45,-9.2>
+v*s:<4.5,8.25,300.9>
+s*v:<4.5,8.25,300.9>
+?*=?:<4.5,8.25,300.9>
+?*?(vec):<2.,6.,13.5>
+?*=?(vec):<2.,6.,13.5>
+?/?:<0.666667,-0.0333333,15.>
+?/=?:<0.666667,-0.0333333,15.>
+?/?(vec):<1.,-0.0333333,0.666667>
+?/=?(vec):<1.,-0.0333333,0.666667>
+dot_1:0.
+dot_2:37.8
+cross:<-15.,-2.,39.>
+length:3.74166
+length_squared:14.
+distance:7.48599
+normalize:<0.872964,-0.465581,0.145494>
+project:<0.,0.,0.>
+project_2:<5.88034,0.588034,2.35214>
+reflect:<5.28917,6.,-1.57084>
+refract:<0.613403,0.736083,0.286214>
+refract:<0.461205,0.82471,-0.327329>
+faceforward_nochange:<4.,5.5,2.>
+faceforward_flip:<-4.,-5.5,-2.>
Index: tests/vector_math/.expect/vec3_int.txt
===================================================================
--- tests/vector_math/.expect/vec3_int.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/.expect/vec3_int.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,37 @@
+ctor(x,y):<1,2,3>
+copy ctor:<1,2,3>
+?/?(vec)<5,-1,1>
+?/=?(vec)<5,-1,1>
+++?<2,-2,101>
+--?<0,-4,99>
+?++(part 1)<1,-3,100>
+?++(part 2)<2,-2,101>
+?--(part 1)<1,-3,100>
+?--(part 2)<0,-4,99>
+?%?(int)<0,1,0>
+?%=?(int)<0,1,0>
+?%?(vec)<-2,3,9>
+?%=?(vec)<-2,3,9>
+?&?(int)<0,0,8>
+?&=?(int)<0,0,8>
+?&?(vec)<1,6,1>
+?&=?(vec)<1,6,1>
+?|?(int)<7,3,7>
+?|=?(int)<7,3,7>
+?|?(vec)<7,1,7>
+?|=?(vec)<7,1,7>
+?^?(int)<12,0,13>
+?^=?(int)<12,0,13>
+?^?(vec)<12,5,9>
+?^=?(vec)<12,5,9>
+?<<?(int)<4,2,8>
+?<<=?(int)<4,2,8>
+?<<?(vec)<32,4,2>
+?<<=?(vec)<32,4,2>
+?>>?(int)<1,0,2>
+?>>=?(int)<1,0,2>
+?>>?(vec)<2,1,2>
+?>>=?(vec)<2,1,2>
+~?<-6,-3,-2>
+dot:4
+length_squared:29
Index: tests/vector_math/.expect/vec4_float.txt
===================================================================
--- tests/vector_math/.expect/vec4_float.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/.expect/vec4_float.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,30 @@
+ctor(x,y):<1.,2.,3.,4.>
+copy ctor:<1.,2.,3.,4.>
+assignment:<3.,4.2,-2.,-100.>
+move assignment:<1.,2.,3.,4.>
+zero-init:<0.,0.,0.,0.>
+zero-assign:<0.,0.,0.,0.>
+fill-ctor:<1.23,1.23,1.23,1.23>
+?-?:<0.02,0.43,-0.999998,-1e-06.>
+?-=?:<0.02,0.43,-0.999998,-1e-06.>
+-?:<-0.02,-0.43,0.999998,1e-06.>
+?+?:<2.3,2.45,-9.2,-12.5>
+?+=?:<2.3,2.45,-9.2,-12.5>
+v*s:<4.5,8.25,300.9,-6.>
+s*v:<4.5,8.25,300.9,-6.>
+?*=?:<4.5,8.25,300.9,-6.>
+?/?:<0.666667,-0.0333333,15.,-5.>
+?/=?:<0.666667,-0.0333333,15.,-5.>
+dot_1:0.
+dot_2:47.7345
+length:5.47723
+length_squared:30.
+distance:10.9563
+normalize:<0.433872,-0.231399,0.072312,0.867745>
+project:<-0.,-0.,-0.,-1.>
+project_2:<2.91736,0.291736,1.16694,0.583471>
+reflect:<5.26163,6.,-1.19906,1.52327>
+refract:<0.608838,0.730606,0.284084,0.121768>
+refract:<0.457773,0.822893,-0.333534,-0.0452281>
+faceforward_nochange:<4.,5.5,2.,2.>
+faceforward_flip:<4.,5.5,2.,2.>
Index: tests/vector_math/.expect/vec4_int.txt
===================================================================
--- tests/vector_math/.expect/vec4_int.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/.expect/vec4_int.txt	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,37 @@
+ctor(x,y):<1,2,3,4>
+copy ctor:<1,2,3,4>
+?/?(vec)<5,-1,1,5>
+?/=?(vec)<5,-1,1,5>
+++?<2,-2,101,2>
+--?<0,-4,99,2>
+?++(part 1)<1,-3,100,1>
+?++(part 2)<2,-2,101,2>
+?--(part 1)<1,-3,100,3>
+?--(part 2)<0,-4,99,2>
+?%?(int)<0,1,0,1>
+?%=?(int)<0,1,0,1>
+?%?(vec)<-2,3,9,4>
+?%=?(vec)<-2,3,9,4>
+?&?(int)<0,0,8,8>
+?&=?(int)<0,0,8,8>
+?&?(vec)<1,6,1,5>
+?&=?(vec)<1,6,1,5>
+?|?(int)<7,3,7,10>
+?|=?(int)<7,3,7,10>
+?|?(vec)<7,1,7,10>
+?|=?(vec)<7,1,7,10>
+?^?(int)<12,0,13,8>
+?^=?(int)<12,0,13,8>
+?^?(vec)<12,5,9,10>
+?^=?(vec)<12,5,9,10>
+?<<?(int)<4,2,8,12>
+?<<=?(int)<4,2,8,12>
+?<<?(vec)<32,4,2,8>
+?<<=?(vec)<32,4,2,8>
+?>>?(int)<1,0,2,3>
+?>>=?(int)<1,0,2,3>
+?>>?(vec)<2,1,2,4>
+?>>=?(vec)<2,1,2,4>
+~?<-6,-3,-2,0>
+dot:-11
+length_squared:54
Index: tests/vector_math/glm_equivalents/vec2_float.cc
===================================================================
--- tests/vector_math/glm_equivalents/vec2_float.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/glm_equivalents/vec2_float.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,28 @@
+#include <glm/glm.hpp>
+#include <iostream>
+
+using namespace glm;
+
+std::ostream& operator<<(std::ostream& os, glm::tvec2<float>& v) {
+    os << "<" << v.x << "," << v.y << ">";
+    return os;
+}
+std::ostream& operator<<(std::ostream& os, glm::tvec2<float>&& v) {
+    os << "<" << v.x << "," << v.y << ">";
+    return os;
+}
+
+float length_squared(glm::tvec2<float> v) {
+    return glm::length(v) * glm::length(v);
+}
+
+tvec2<float> project(glm::tvec2<float> u, glm::tvec2<float> v) {
+    return normalize(v) * dot(u, normalize(v));
+}
+
+int main() {
+    tvec2<float> v1 = {5.f,6.f};
+    tvec2<float> v2 = {0.f,-1.f};
+    std::cout << "refract:" << refract(v1,normalize(v2),1.f) << std::endl;
+    std::cout << "refract:" << refract(v1,normalize(v2),1.f/1.33f) << std::endl;
+}
Index: tests/vector_math/glm_equivalents/vec2_int.cc
===================================================================
--- tests/vector_math/glm_equivalents/vec2_int.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/glm_equivalents/vec2_int.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,120 @@
+#include <glm/glm.hpp>
+#include <iostream>
+
+using namespace glm;
+
+std::ostream& operator<<(std::ostream& os, glm::tvec2<int>& v) {
+    os << "<" << v.x << "," << v.y << ">";
+    return os;
+}
+std::ostream& operator<<(std::ostream& os, glm::tvec2<int>&& v) {
+    os << "<" << v.x << "," << v.y << ">";
+    return os;
+}
+
+int length_squared(glm::tvec2<float> v) {
+    return dot(v, v);
+}
+
+// tvec2<int> project(glm::tvec2<int> u, glm::tvec2<int> v) {
+//     return normalize(v) * dot(u, normalize(v));
+// }
+
+int main(void) {
+    tvec2<int> v1 = {1,2};
+    std::cout << "ctor(x,y):" << v1 << std::endl;
+
+    tvec2<int> v2 = v1;
+    std::cout << "copy ctor:" << v2 << std::endl;
+
+    v2 = (tvec2<int>){3, 4};
+    std::cout << "assignment:" << v2 << std::endl;
+
+    v2 = v1;
+    std::cout << "move assignment:" << v2 << std::endl;
+
+
+                // NOT SUPPORTED BY GLM
+    // tvec2<int> v3 = 0;
+    // std::cout << "zero-init:" << v3 << std::endl;
+
+    // v1 = 0;
+    // std::cout << "zero-assign:" << v1 << std::endl;
+
+    // tvec2<int> v4 = {1.23f};
+    // std::cout << "fill-ctor:" << v4 << std::endl;
+
+
+    v1 = (tvec2<int>){1, 3};
+    std::cout << "?-?:" << (v1 - (tvec2<int>){1003, -14}) << std::endl;
+
+    v1 -= (tvec2<int>){1003, -14};
+    std::cout << "?-=?:" << v1 << std::endl;
+
+    v1 = -v1;
+    std::cout << "-?:" << v1 << std::endl;
+
+    v1 = (tvec2<int>){1, 3};
+    std::cout << "?+?:" << (v1 + (tvec2<int>){1003, -14}) << std::endl;
+
+    v1 += (tvec2<int>){1003, -14};
+    std::cout << "?+=?:" << v1 << std::endl;
+
+    v1 = (tvec2<int>){15, 275};
+    std::cout << "v*s:" << v1 * 3 << std::endl;
+
+    std::cout << "s*v:" << 3 * v1 << std::endl;
+
+    v1 *= 3;
+    std::cout << "?*=?:" << v1 << std::endl;
+
+    v1 = (tvec2<int>){21, -10};
+    std::cout << "?/?:" << (v1 / 3) << std::endl;
+
+    v1 /= 3;
+    std::cout << "?/=?:" << v1 << std::endl;
+
+    //      FORCE GLM TO FLOAT, BUT WORKS
+    tvec2<float> v1_f = (tvec2<float>){2.f, 3};
+    tvec2<float> v2_f = (tvec2<float>){-3.f, 2.f};
+    std::cout << "dot_1:" << dot(v1_f, v2_f) << std::endl;
+
+    v2_f = (tvec2<float>){13, 2};
+    std::cout << "dot_2:" << dot(v1_f, v2_f) << std::endl;
+
+
+    //      NOT IN GLM BUT
+    //      USES DOT, BUT THAT'S FLOATING-POINT ONLY!
+    std::cout << "length_squared:" << length_squared(v1_f) << std::endl;
+
+    //      FLOATING-POINT ONLY GLM
+    // v1 = (tvec2<int>){100, -101};
+    // v2 = (tvec2<int>){6, 3};
+    // std::cout << "distance:" << distance(v1, v2) << std::endl;
+    //
+    // std::cout << "normalize:" << normalize(v2) << std::endl;
+    //
+    // v1 = (tvec2<int>){1,0};
+    // std::cout << "normalize_2:" << normalize(v1) << std::endl;
+    //
+    // std::cout << "project:" << project((tvec2<int>){5,6}, (tvec2<int>){1, 0}) << std::endl;
+    // std::cout << "project_2:" << project((tvec2<int>){5,6}, (tvec2<int>){1, 1}) << std::endl;
+    //
+    // v1 = (tvec2<int>){5,6};
+    // v2 = (tvec2<int>){1,0};
+    // std::cout << "reflect:" << reflect(v1,v2) << std::endl;
+    //
+    // v2 = (tvec2<int>){0,-1};
+    // std::cout << "refract:" << refract(v1,normalize(v2),1) << std::endl;
+    // std::cout << "refract:" << refract(v1,normalize(v2),1.f/1.33f) << std::endl;
+    //
+    //
+    //      FLOATING-POINT ONLY IN GLM BUT JUST USES DOT SO...
+    tvec2<float> geometric_normal = {5,6};
+    tvec2<float> perturbed_normal = {4,5};
+    tvec2<float> eyeline = {-1,0};
+    std::cout << "faceforward_nochange:" << faceforward(perturbed_normal, eyeline, geometric_normal) << std::endl;
+
+    eyeline = (tvec2<int>){1,0};
+    std::cout << "faceforward_flip:" << faceforward(perturbed_normal, eyeline, geometric_normal) << std::endl;
+}
Index: tests/vector_math/glm_equivalents/vec3_float.cc
===================================================================
--- tests/vector_math/glm_equivalents/vec3_float.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/glm_equivalents/vec3_float.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,117 @@
+#include <glm/glm.hpp>
+#include <iostream>
+#include <ostream>
+
+using namespace glm;
+
+std::ostream& operator<<(std::ostream& os, glm::tvec3<float>& v) {
+    os << "<" << v.x << "," << v.y << "," << v.z << ">";
+    return os;
+}
+std::ostream& operator<<(std::ostream& os, glm::tvec3<float>&& v) {
+    os << "<" << v.x << "," << v.y << "," << v.z << ">";
+    return os;
+}
+
+float length_squared(glm::tvec3<float> v) {
+    return glm::length(v) * glm::length(v);
+}
+
+tvec3<float> project(glm::tvec3<float> u, glm::tvec3<float> v) {
+    return normalize(v) * dot(u, normalize(v));
+}
+
+int main(void) {
+    tvec3<float> v1 = {1.f,2.f,3.f};
+    std::cout << "ctor(x,y):" << v1 << std::endl;
+
+    tvec3<float> v2 = v1;
+    std::cout << "copy ctor:" << v2 << std::endl;
+
+    v2 = (tvec3<float>){3.f, 4.2f, -2.f};
+    std::cout << "assignment:" << v2 << std::endl;
+
+    v2 = v1;
+    std::cout << "move assignment:" << v2 << std::endl;
+
+    // tvec3<float> v3 = 0;
+    // std::cout << "zero-init:" << v3 << std::endl;
+
+    // v1 = 0;
+    // std::cout << "zero-assign:" << v1 << std::endl;
+    //
+    // tvec3<float> v4 = {1.23f};
+    // std::cout << "fill-ctor:" << v4 << std::endl;
+    //
+    v1 = (tvec3<float>){1.23f, 3.43f, 0.000002f};
+    std::cout << "?-?:" << (v1 - (tvec3<float>){1.21f,3,1}) << std::endl;
+
+    v1 -= (tvec3<float>){1.23f, 3.43f, 0.000002f};
+    std::cout << "?-=?:" << v1 << std::endl;
+
+    v1 = -v1;
+    std::cout << "-?:" << v1 << std::endl;
+
+
+
+
+
+    v1 = (tvec3<float>){1.5f, 2.75f, -14.2f};
+    std::cout << "?+?:" << (v1 + (tvec3<float>){0.8f, -0.3f, 5}) << std::endl;
+
+    v1 += (tvec3<float>){0.8f, -0.3f, 5};
+    std::cout << "?+=?:" << v1 << std::endl;
+
+    v1 = (tvec3<float>){1.5f, 2.75f, 100.3f};
+    std::cout << "v*s:" << v1 * 3.f << std::endl;
+
+    std::cout << "s*v:" << 3.f * v1 << std::endl;
+
+    v1 *= 3;
+    std::cout << "?*=?:" << v1 << std::endl;
+
+    v1 = (tvec3<float>){2, -0.1f, 45};
+    std::cout << "?/?:" << (v1 / 3.f) << std::endl;
+
+    v1 /= 3.f;
+    std::cout << "?/=?:" << v1 << std::endl;
+
+
+    v1 = (tvec3<float>){4,2,3};
+    v2 = (tvec3<float>){0,-3,2};
+    std::cout << "dot_1:" << dot(v1, v2) << std::endl;
+
+    v2 = (tvec3<float>){1.3f, -2, 12.2};
+    std::cout << "dot_2:" << dot(v1, v2) << std::endl;
+
+
+    v1 = (tvec3<float>){1,2,3};
+    std::cout << "length:" << length(v1) << std::endl;
+
+    std::cout << "length_squared:" << length_squared(v1) << std::endl;
+
+    v2 = (tvec3<float>){6, -3.2f, 1};
+    std::cout << "distance:" << distance(v1, v2) << std::endl;
+
+    std::cout << "normalize:" << normalize(v2) << std::endl;
+
+    std::cout << "project:" << project((tvec3<float>){5,6,0}, (tvec3<float>){0, 0, 1}) << std::endl;
+    std::cout << "project_2:" << project((tvec3<float>){5,6,3.2f}, (tvec3<float>){10, 1, 4}) << std::endl;
+
+    v1 = (tvec3<float>){5,6,2.333f};
+    v2 = (tvec3<float>){1,0,-13.5f};
+    std::cout << "reflect:" << reflect(v1,normalize(v2)) << std::endl;
+
+    v2 = (tvec3<float>){0,-1,2};
+    std::cout << "refract:" << refract(normalize(v1),normalize(v2),1.f) << std::endl;
+    std::cout << "refract:" << refract(normalize(v1),normalize(v2),1.f/1.33f) << std::endl;
+
+    tvec3<float> geometric_normal = {5,6,1};
+    tvec3<float> perturbed_normal = {4,5.5f,2};
+    tvec3<float> eyeline = {-1,0.002f,-1.0345f};
+    std::cout << "faceforward_nochange:" << faceforward(perturbed_normal, eyeline, geometric_normal) << std::endl;
+
+    eyeline = (tvec3<float>){1,0.002f,-1.0345f};
+    std::cout << "faceforward_flip:" << faceforward(perturbed_normal, eyeline, geometric_normal) << std::endl;
+}
+
Index: tests/vector_math/glm_equivalents/vec4_float.cc
===================================================================
--- tests/vector_math/glm_equivalents/vec4_float.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/glm_equivalents/vec4_float.cc	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,114 @@
+#include <glm/glm.hpp>
+#include <iostream>
+#include <ostream>
+
+using namespace glm;
+
+std::ostream& operator<<(std::ostream& os, glm::tvec4<float>& v) {
+    os << "<" << v.x << "," << v.y << "," << v.z << "," << v.w << ">";
+    return os;
+}
+std::ostream& operator<<(std::ostream& os, glm::tvec4<float>&& v) {
+    os << "<" << v.x << "," << v.y << "," << v.z << "," << v.w << ">";
+    return os;
+}
+
+float length_squared(glm::tvec4<float> v) {
+    return glm::length(v) * glm::length(v);
+}
+
+tvec4<float> project(glm::tvec4<float> u, glm::tvec4<float> v) {
+    return normalize(v) * dot(u, normalize(v));
+}
+
+int main(void) {
+    tvec4<float> v1 = {1.f,2.f,3.f,4.f};
+    std::cout << "ctor(x,y):" << v1 << std::endl;
+
+    tvec4<float> v2 = v1;
+    std::cout << "copy ctor:" << v2 << std::endl;
+
+    v2 = (tvec4<float>){3.f,4.2f,-2.f,-100.f};
+    std::cout << "assignment:" << v2 << std::endl;
+
+    v2 = v1;
+    std::cout << "move assignment:" << v2 << std::endl;
+
+    // tvec4<float> v3 = 0;
+    // std::cout << "zero-init:" << v3 << std::endl;
+
+    // v1 = 0;
+    // std::cout << "zero-assign:" << v1 << std::endl;
+    //
+    // tvec4<float> v4 = {1.23f};
+    // std::cout << "fill-ctor:" << v4 << std::endl;
+    //
+    v1 = (tvec4<float>){1.23f, 3.43f, 0.000002f, -0.000002f};
+    std::cout << "?-?:" << (v1 - (tvec4<float>){1.21f,3,1,-0.000001f}) << std::endl;
+
+    v1 -= (tvec4<float>){1.21f,3,1,-0.000001f};
+    std::cout << "?-=?:" << v1 << std::endl;
+
+    v1 = -v1;
+    std::cout << "-?:" << v1 << std::endl;
+
+
+    v1 = (tvec4<float>){1.5f, 2.75f, -14.2f, -13.5f};
+    std::cout << "?+?:" << (v1 + (tvec4<float>){0.8f, -0.3f, 5, 1}) << std::endl;
+
+    v1 += (tvec4<float>){0.8f, -0.3f, 5, 1};
+    std::cout << "?+=?:" << v1 << std::endl;
+
+    v1 = (tvec4<float>){1.5f, 2.75f, 100.3f, -2.f};
+    std::cout << "v*s:" << v1 * 3.f << std::endl;
+
+    std::cout << "s*v:" << 3.f * v1 << std::endl;
+
+    v1 *= 3;
+    std::cout << "?*=?:" << v1 << std::endl;
+
+    v1 = (tvec4<float>){2, -0.1f, 45, -15};
+    std::cout << "?/?:" << (v1 / 3.f) << std::endl;
+
+    v1 /= 3.f;
+    std::cout << "?/=?:" << v1 << std::endl;
+
+
+    v1 = (tvec4<float>){4,3,2,1};
+    v2 = (tvec4<float>){0,4,-3,-6};
+    std::cout << "dot_1:" << dot(v1, v2) << std::endl;
+
+    v2 = (tvec4<float>){1.3f, -2, 12.2, 24.1345};
+    std::cout << "dot_2:" << dot(v1, v2) << std::endl;
+
+
+    v1 = (tvec4<float>){1,2,3,4};
+    std::cout << "length:" << length(v1) << std::endl;
+
+    std::cout << "length_squared:" << length_squared(v1) << std::endl;
+
+    v2 = (tvec4<float>){6, -3.2f, 1, 12};
+    std::cout << "distance:" << distance(v1, v2) << std::endl;
+
+    std::cout << "normalize:" << normalize(v2) << std::endl;
+
+    std::cout << "project:" << project((tvec4<float>){5,6,0,-1}, (tvec4<float>){0, 0, 0, 1}) << std::endl;
+    std::cout << "project_2:" << project((tvec4<float>){5,6,3.2f,-16.75f}, (tvec4<float>){10, 1, 4, 2}) << std::endl;
+
+    v1 = (tvec4<float>){5,6,2.333f,1};
+    v2 = (tvec4<float>){1,0,-13.5f,2};
+    std::cout << "reflect:" << reflect(v1,normalize(v2)) << std::endl;
+
+    v2 = (tvec4<float>){0,-1,2,0.5f};
+    std::cout << "refract:" << refract(normalize(v1),normalize(v2),1.f) << std::endl;
+    std::cout << "refract:" << refract(normalize(v1),normalize(v2),1.f/1.33f) << std::endl;
+
+    tvec4<float> geometric_normal = {5,6,1,2};
+    tvec4<float> perturbed_normal = {4,5.5f,2,2};
+    tvec4<float> eyeline = {-1,0.002f,-1.0345f,-2};
+    std::cout << "faceforward_nochange:" << faceforward(perturbed_normal, eyeline, geometric_normal) << std::endl;
+
+    eyeline = (tvec4<float>){1,0.002f,-1.0345f,-2};
+    std::cout << "faceforward_flip:" << faceforward(perturbed_normal, eyeline, geometric_normal) << std::endl;
+}
+
Index: tests/vector_math/vec2_double.cfa
===================================================================
--- tests/vector_math/vec2_double.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/vec2_double.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,12 @@
+#include <vec/vec2.hfa>
+#include <fstream.hfa>
+#include <limits.hfa>
+
+int main(void) {
+    vec2(double) v1 = {(double)MAX,(double)MIN};
+    sout | "ctor(x,y):" | v1;
+
+    vec2(double) v2 = {3.4144423, -1.2342315};
+    vec2(double) v3 = {100.23125123, 14.23};
+    sout | "dot(v2,v3):" | dot(v2, v3);
+}
Index: tests/vector_math/vec2_float.cfa
===================================================================
--- tests/vector_math/vec2_float.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/vec2_float.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,127 @@
+#include <vec/vec2.hfa>
+#include <fstream.hfa>
+
+int main(void) {
+    vec2(float) v1 = {1.f,2.f};
+    sout | "ctor(x,y):" | v1;
+
+    vec2(float) v2 = v1;
+    sout | "copy ctor:" | v2;
+
+    v2 = (vec2(float)){3.f, 4.2f};
+    sout | "assignment:" | v2;
+
+    v2 = v1;
+    sout | "move assignment:" | v2;
+
+    vec2(float) v3 = 0;
+    sout | "zero-init:" | v3;
+
+    v1 = 0;
+    sout | "zero-assign:" | v1;
+
+    vec2(float) v4 = {1.23f};
+    sout | "fill-ctor:" | v4;
+
+    v1 = (vec2(float)){1.23f, 3.43f};
+    sout | "?-?:" | (v1 - (vec2(float)){1.21f,3});
+
+    v1 -= (vec2(float)){-2, 10};
+    sout | "?-=?:" | v1;
+
+    v1 = -v1;
+    sout | "-?:" | v1;
+
+    v1 = (vec2(float)){1.5f, 2.75f};
+    sout | "?+?:" | (v1 + (vec2(float)){0.8f, -0.3f});
+
+    v1 += (vec2(float)){0.8f, -0.3f};
+    sout | "?+=?:" | v1;
+
+    v1 = (vec2(float)){1.5f, 2.75f};
+    sout | "v*s:" | v1 * 3.f;
+
+    sout | "s*v:" | 3.f * v1;
+
+    v1 *= 3;
+    sout | "?*=?:" | v1;
+
+    {
+        vec2(float) u = {1.2f, -5.8f};
+        vec2(float) v = {2.2f, 3};
+        sout | "?*?(vec)" | (u * v);
+        sout | "?*=?(vec)" | (u *= v);
+    }
+
+
+    v1 = (vec2(float)){2, -0.1f};
+    sout | "?/?:" | (v1 / 3);
+
+    v1 /= 3;
+    sout | "?/=?:" | v1;
+
+    {
+        vec2(float) u = {1.2f, -5.8f};
+        vec2(float) v = {2.2f, 3};
+        sout | "?/?(vec)" | (u / v);
+        sout | "?/=?(vec)" | (u /= v);
+    }
+
+    {
+        vec2(float) u = {1.2f, -3.5f};
+        sout | "++?" | ++u;
+    }
+    {
+        vec2(float) u = {1.2f, -3.5f};
+        sout | "--?" | ++u;
+    }
+    {
+        vec2(float) u = {1.2f, -3.5f};
+        sout | "?++(part 1)" | u++;
+        sout | "?++(part 2)" | u;
+    }
+    {
+        vec2(float) u = {1.2f, -3.5f};
+        sout | "?--(part 1)" | u--;
+        sout | "?--(part 2)" | u;
+    }
+
+    v1 = (vec2(float)){2, 3};
+    v2 = (vec2(float)){-3, 2};
+    sout | "dot_1:" | dot(v1, v2);
+
+    v2 = (vec2(float)){1.3f, 2};
+    sout | "dot_2:" | dot(v1, v2);
+
+    v1 = (vec2(float)){4, 3};
+    sout | "length:" | length(v1);
+
+    sout | "length_squared:" | length_squared(v1);
+
+    v2 = (vec2(float)){6, -3.2f};
+    sout | "distance:" | distance(v1, v2);
+
+    sout | "normalize:" | normalize(v2);
+
+    v1 = (vec2(float)){1,0};
+    sout | "normalize_2:" | normalize(v1);
+
+    sout | "project:" | project((vec2(float)){5,6}, (vec2(float)){1, 0});
+    sout | "project_2:" | project((vec2(float)){5,6}, (vec2(float)){1, 1});
+
+    v1 = (vec2(float)){5,6};
+    v2 = (vec2(float)){1,0};
+    sout | "reflect:" | reflect(v1,v2);
+
+    v2 = (vec2(float)){0,-1};
+    sout | "refract:" | refract(v1,normalize(v2),1.f);
+    sout | "refract:" | refract(v1,normalize(v2),1.f/1.33f);
+
+    vec2(float) geometric_normal = {5,6};
+    vec2(float) perturbed_normal = {4,5.5f};
+    vec2(float) eyeline = {-1,0.002f};
+    sout | "faceforward_nochange:" | faceforward(perturbed_normal, eyeline, geometric_normal);
+
+    eyeline = (vec2(float)){1,0.002f};
+    sout | "faceforward_flip:" | faceforward(perturbed_normal, eyeline, geometric_normal);
+}
Index: tests/vector_math/vec2_int.cfa
===================================================================
--- tests/vector_math/vec2_int.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/vec2_int.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,191 @@
+#include <vec/vec2.hfa>
+#include <fstream.hfa>
+#include <limits.hfa>
+
+int main(void) {
+    vec2(int) v1 = {1,2};
+    sout | "ctor(x,y):" | v1;
+
+    vec2(int) v2 = v1;
+    sout | "copy ctor:" | v2;
+
+    v2 = (vec2(int)){3, 4};
+    sout | "assignment:" | v2;
+
+    v2 = v1;
+    sout | "move assignment:" | v2;
+
+    vec2(int) v3 = 0;
+    sout | "zero-init:" | v3;
+
+    v1 = 0;
+    sout | "zero-assign:" | v1;
+
+    vec2(int) v4 = {123};
+    sout | "fill-ctor:" | v4;
+
+    v1 = (vec2(int)){1, 3};
+    sout | "?-?:" | (v1 - (vec2(int)){1003, -14});
+
+    v1 -= (vec2(int)){1003, -14};
+    sout | "?-=?:" | v1;
+
+    v1 = -v1;
+    sout | "-?:" | v1;
+
+    v1 = (vec2(int)){1, 3};
+    sout | "?+?:" | (v1 + (vec2(int)){1003, -14});
+
+    v1 += (vec2(int)){1003, -14};
+    sout | "?+=?:" | v1;
+
+    v1 = (vec2(int)){15, 275};
+    sout | "v*s:" | v1 * 3;
+
+    sout | "s*v:" | 3 * v1;
+
+    v1 *= 3;
+    sout | "?*=?:" | v1;
+
+    {
+        vec2(int) u = {1, -5};
+        vec2(int) v = {2, 3};
+        sout | "?*?(vec)" | (u * v);
+        sout | "?*=?(vec)" | (u *= v);
+    }
+
+    v1 = (vec2(int)){21, -10};
+    sout | "?/?:" | (v1 / 3);
+
+    v1 /= 3;
+    sout | "?/=?:" | v1;
+
+    {
+        vec2(int) u = {11, -5};
+        vec2(int) v = {2, 3};
+        sout | "?/?(vec)" | (u / v);
+        sout | "?/=?(vec)" | (u /= v);
+    }
+
+    {
+        vec2(int) u = {1, -3};
+        sout | "++?" | ++u;
+    }
+    {
+        vec2(int) u = {1, -3};
+        sout | "--?" | ++u;
+    }
+    {
+        vec2(int) u = {1, -3};
+        sout | "?++(part 1)" | u++;
+        sout | "?++(part 2)" | u;
+    }
+    {
+        vec2(int) u = {1, -3};
+        sout | "?--(part 1)" | u--;
+        sout | "?--(part 2)" | u;
+    }
+
+    // %
+    {
+        vec2(int) u = {-12, 123};
+        sout | "?%?(int)" | (u % 2);
+        sout | "?%=?(int)" | (u %= 2);
+    }
+    {
+        vec2(int) u = {-12, 27};
+        vec2(int) v = {5, 4};
+        sout | "?%?(vec)" | (u % v);
+        sout | "?%=?(vec)" | (u %= v);
+    }
+
+
+    // &
+    {
+        vec2(int) u = {5, 8};
+        sout | "?&?(int)" | (u & (1 << 3));
+        sout | "?&=?(int)" | (u &= (1 << 3));
+    }
+    {
+        vec2(int) u = {5, 6};
+        vec2(int) v = {3, 7};
+        sout | "?&?(vec)" | (u & v);
+        sout | "?&=?(vec)" | (u &= v);
+    }
+
+    // |
+    {
+        vec2(int) u = {5, 1};
+        sout | "?|?(int)" | (u | 2);
+        sout | "?|=?(int)" | (u |= 2);
+    }
+    {
+        vec2(int) u = {5, 1};
+        vec2(int) v = {2, 1};
+        sout | "?|?(vec)" | (u | v);
+        sout | "?|=?(vec)" | (u |= v);
+    }
+
+    // ^
+    {
+        vec2(int) u = {9, 5};
+        sout | "?^?(int)" | (u ^ 5);
+        sout | "?^=?(int)" | (u ^= 5);
+    }
+    {
+        vec2(int) u = {9, 5};
+        vec2(int) v = {5, 0};
+        sout | "?^?(vec)" | (u ^ v);
+        sout | "?^=?(vec)" | (u ^= v);
+    }
+
+    // <<
+    {
+        vec2(int) u = {2, 1};
+        sout | "?<<?(int)" | (u << 1);
+        sout | "?<<=?(int)" | (u <<= 1);
+    }
+    {
+        vec2(int) u = {1, 1};
+        vec2(int) v = {5, 2};
+        sout | "?<<?(vec)" | (u << v);
+        sout | "?<<=?(vec)" | (u <<= v);
+    }
+
+    // >>
+    {
+        vec2(int) u = {2, 1};
+        sout | "?>>?(int)" | (u >> 1);
+        sout | "?>>=?(int)" | (u >>= 1);
+    }
+    {
+        vec2(int) u = {5, 7};
+        vec2(int) v = {1, 2};
+        sout | "?>>?(vec)" | (u >> v);
+        sout | "?>>=?(vec)" | (u >>= v);
+    }
+
+    // ~
+    {
+        vec2(int) u = {5, 2};
+        sout | "~?" | (~u);
+    }
+
+    v1 = (vec2(int)){2, 3};
+    v2 = (vec2(int)){-3, 2};
+    sout | "dot_1:" | dot(v1, v2);
+
+    v2 = (vec2(int)){13, 2};
+    sout | "dot_2:" | dot(v1, v2);
+
+    v1 = (vec2(int)){4, 3};
+    sout | "length_squared:" | length_squared(v1);
+
+    vec2(int) geometric_normal = {5,6};
+    vec2(int) perturbed_normal = {4,5};
+    vec2(int) eyeline = {-1,0};
+    sout | "faceforward_nochange:" | faceforward(perturbed_normal, eyeline, geometric_normal);
+
+    eyeline = (vec2(int)){1,0};
+    sout | "faceforward_flip:" | faceforward(perturbed_normal, eyeline, geometric_normal);
+}
Index: tests/vector_math/vec2_ldouble.cfa
===================================================================
--- tests/vector_math/vec2_ldouble.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/vec2_ldouble.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,12 @@
+#include <vec/vec2.hfa>
+#include <fstream.hfa>
+#include <limits.hfa>
+
+int main(void) {
+    vec2(long double) v1 = {(long double)MAX,(long double)MIN};
+    sout | "ctor(x,y):" | v1;
+
+    vec2(long double) v2 = {3.4144423, -1.2342315};
+    vec2(long double) v3 = {100.23125123, 14.23};
+    sout | "dot(v2,v3):" | dot(v2, v3);
+}
Index: tests/vector_math/vec2_uint.cfa
===================================================================
--- tests/vector_math/vec2_uint.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/vec2_uint.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,12 @@
+#include <vec/vec2.hfa>
+#include <fstream.hfa>
+#include <limits.hfa>
+
+int main(void) {
+    vec2(unsigned) v1 = {(unsigned)MAX,(unsigned)MIN};
+    sout | "ctor(x,y):" | v1;
+
+    vec2(unsigned) v2 = {3, 4};
+    vec2(unsigned) v3 = {100, 13};
+    sout | "dot(v2,v3):" | dot(v2, v3);
+}
Index: tests/vector_math/vec3_float.cfa
===================================================================
--- tests/vector_math/vec3_float.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/vec3_float.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,107 @@
+#include <vec/vec3.hfa>
+#include <fstream.hfa>
+
+int main(void) {
+    vec3(float) v1 = {1.f,2.f,3.f};
+    sout | "ctor(x,y):" | v1;
+
+    vec3(float) v2 = v1;
+    sout | "copy ctor:" | v2;
+
+    v2 = (vec3(float)){3.f, 4.2f, -2.f};
+    sout | "assignment:" | v2;
+
+    v2 = v1;
+    sout | "move assignment:" | v2;
+
+    vec3(float) v3 = 0;
+    sout | "zero-init:" | v3;
+
+    v1 = 0;
+    sout | "zero-assign:" | v1;
+
+    vec3(float) v4 = {1.23f};
+    sout | "fill-ctor:" | v4;
+
+    v1 = (vec3(float)){1.23f, 3.43f, 0.000002f};
+    sout | "?-?:" | (v1 - (vec3(float)){1.21f,3,1});
+
+    v1 -= (vec3(float)){1.23f, 3.43f, 0.000002f};
+    sout | "?-=?:" | v1;
+
+    v1 = -v1;
+    sout | "-?:" | v1;
+
+    v1 = (vec3(float)){1.5f, 2.75f, -14.2f};
+    sout | "?+?:" | (v1 + (vec3(float)){0.8f, -0.3f, 5});
+
+    v1 += (vec3(float)){0.8f, -0.3f, 5};
+    sout | "?+=?:" | v1;
+
+    v1 = (vec3(float)){1.5f, 2.75f, 100.3f};
+    sout | "v*s:" | v1 * 3.f;
+
+    sout | "s*v:" | 3.f * v1;
+
+    v1 *= 3;
+    sout | "?*=?:" | v1;
+
+    {
+        vec3(float) u = {1, 2, 3};
+        vec3(float) v = {2, 3, 4.5f};
+        sout | "?*?(vec):" | (u * v);
+        sout | "?*=?(vec):" | (u *= v);
+    }
+
+    v1 = (vec3(float)){2, -0.1f, 45};
+    sout | "?/?:" | (v1 / 3.f);
+
+    v1 /= 3.f;
+    sout | "?/=?:" | v1;
+
+    {
+        vec3(float) u = {2, -0.1f, 3};
+        vec3(float) v = {2, 3, 4.5f};
+        sout | "?/?(vec):" | (u / v);
+        sout | "?/=?(vec):" | (u /= v);
+    }
+
+    v1 = (vec3(float)){4,2,3};
+    v2 = (vec3(float)){0,-3,2};
+    sout | "dot_1:" | dot(v1, v2);
+
+    v2 = (vec3(float)){1.3f, -2, 12.2f};
+    sout | "dot_2:" | dot(v1, v2);
+
+    v2 = cross((vec3(float)){3,-3,1},(vec3(float)){4,9,2});
+    sout | "cross:" | v2;
+
+    v1 = (vec3(float)){1,2,3};
+    sout | "length:" | length(v1);
+
+    sout | "length_squared:" | length_squared(v1);
+
+    v2 = (vec3(float)){6, -3.2f, 1};
+    sout | "distance:" | distance(v1, v2);
+
+    sout | "normalize:" | normalize(v2);
+
+    sout | "project:" | project((vec3(float)){5,6,0}, (vec3(float)){0, 0, 1});
+    sout | "project_2:" | project((vec3(float)){5,6,3.2f}, (vec3(float)){10, 1, 4});
+
+    v1 = (vec3(float)){5,6,2.333f};
+    v2 = (vec3(float)){1,0,-13.5f};
+    sout | "reflect:" | reflect(v1,v2);
+
+    v2 = (vec3(float)){0,-1,2};
+    sout | "refract:" | refract(normalize(v1),normalize(v2),1.f);
+    sout | "refract:" | refract(normalize(v1),normalize(v2),1.f/1.33f);
+
+    vec3(float) geometric_normal = {5,6,1};
+    vec3(float) perturbed_normal = {4,5.5f,2};
+    vec3(float) eyeline = {-1,0.002f,-1.0345f};
+    sout | "faceforward_nochange:" | faceforward(perturbed_normal, eyeline, geometric_normal);
+
+    eyeline = (vec3(float)){1,0.002f,-1.0345f};
+    sout | "faceforward_flip:" | faceforward(perturbed_normal, eyeline, geometric_normal);
+}
Index: tests/vector_math/vec3_int.cfa
===================================================================
--- tests/vector_math/vec3_int.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/vec3_int.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,130 @@
+#include <vec/vec3.hfa>
+#include <fstream.hfa>
+#include <limits.hfa>
+
+int main(void) {
+    vec3(int) v1 = {1,2,3};
+    sout | "ctor(x,y):" | v1;
+
+    vec3(int) v2 = v1;
+    sout | "copy ctor:" | v2;
+
+    {
+        vec3(int) u = {11, -5, 7};
+        vec3(int) v = {2, 3, 7};
+        sout | "?/?(vec)" | (u / v);
+        sout | "?/=?(vec)" | (u /= v);
+    }
+
+    {
+        vec3(int) u = {1, -3, 100};
+        sout | "++?" | ++u;
+    }
+
+    {
+        vec3(int) u = {1, -3, 100};
+        sout | "--?" | --u;
+    }
+    {
+        vec3(int) u = {1, -3, 100};
+        sout | "?++(part 1)" | u++;
+        sout | "?++(part 2)" | u;
+    }
+    {
+        vec3(int) u = {1, -3, 100};
+        sout | "?--(part 1)" | u--;
+        sout | "?--(part 2)" | u;
+    }
+
+    // %
+    {
+        vec3(int) u = {-12, 123, 100};
+        sout | "?%?(int)" | (u % 2);
+        sout | "?%=?(int)" | (u %= 2);
+    }
+    {
+        vec3(int) u = {-12, 27, 100};
+        vec3(int) v = {5, 4, 13};
+        sout | "?%?(vec)" | (u % v);
+        sout | "?%=?(vec)" | (u %= v);
+    }
+
+
+    // &
+    {
+        vec3(int) u = {5, 1, 8};
+        sout | "?&?(int)" | (u & (1 << 3));
+        sout | "?&=?(int)" | (u &= (1 << 3));
+    }
+    {
+        vec3(int) u = {5, 6, 5};
+        vec3(int) v = {3, 7, 3};
+        sout | "?&?(vec)" | (u & v);
+        sout | "?&=?(vec)" | (u &= v);
+    }
+
+    // |
+    {
+        vec3(int) u = {5, 1, 5};
+        sout | "?|?(int)" | (u | 2);
+        sout | "?|=?(int)" | (u |= 2);
+    }
+    {
+        vec3(int) u = {5, 1, 5};
+        vec3(int) v = {2, 1, 2};
+        sout | "?|?(vec)" | (u | v);
+        sout | "?|=?(vec)" | (u |= v);
+    }
+
+    // ^
+    {
+        vec3(int) u = {9, 5, 8};
+        sout | "?^?(int)" | (u ^ 5);
+        sout | "?^=?(int)" | (u ^= 5);
+    }
+    {
+        vec3(int) u = {9, 5, 8};
+        vec3(int) v = {5, 0, 1};
+        sout | "?^?(vec)" | (u ^ v);
+        sout | "?^=?(vec)" | (u ^= v);
+    }
+
+    // <<
+    {
+        vec3(int) u = {2, 1, 4};
+        sout | "?<<?(int)" | (u << 1);
+        sout | "?<<=?(int)" | (u <<= 1);
+    }
+    {
+        vec3(int) u = {1, 1, 1};
+        vec3(int) v = {5, 2, 1};
+        sout | "?<<?(vec)" | (u << v);
+        sout | "?<<=?(vec)" | (u <<= v);
+    }
+
+    // >>
+    {
+        vec3(int) u = {2, 1, 4};
+        sout | "?>>?(int)" | (u >> 1);
+        sout | "?>>=?(int)" | (u >>= 1);
+    }
+    {
+        vec3(int) u = {5, 7, 8};
+        vec3(int) v = {1, 2, 2};
+        sout | "?>>?(vec)" | (u >> v);
+        sout | "?>>=?(vec)" | (u >>= v);
+    }
+
+    // ~
+    {
+        vec3(int) u = {5, 2, 1};
+        sout | "~?" | (~u);
+    }
+
+    v1 = (vec3(int)){2, 3, 4};
+    v2 = (vec3(int)){13, 2, -7};
+    sout | "dot:" | dot(v1, v2);
+
+    v1 = (vec3(int)){2,3,4};
+    sout | "length_squared:" | length_squared(v1);
+}
Index: tests/vector_math/vec4_float.cfa
===================================================================
--- tests/vector_math/vec4_float.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/vec4_float.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,90 @@
+#include <vec/vec4.hfa>
+#include <fstream.hfa>
+
+int main(void) {
+    vec4(float) v1 = {1.f,2.f,3.f,4.f};
+    sout | "ctor(x,y):" | v1;
+
+    vec4(float) v2 = v1;
+    sout | "copy ctor:" | v2;
+
+    v2 = (vec4(float)){3.f,4.2f,-2.f,-100.f};
+    sout | "assignment:" | v2;
+
+    v2 = v1;
+    sout | "move assignment:" | v2;
+
+    vec4(float) v3 = 0;
+    sout | "zero-init:" | v3;
+
+    v1 = 0;
+    sout | "zero-assign:" | v1;
+
+    vec4(float) v4 = {1.23f};
+    sout | "fill-ctor:" | v4;
+
+    v1 = (vec4(float)){1.23f, 3.43f, 0.000002f, -0.000002f};
+    sout | "?-?:" | (v1 - (vec4(float)){1.21f,3,1,-0.000001f});
+
+    v1 -= (vec4(float)){1.21f,3,1,-0.000001f};
+    sout | "?-=?:" | v1;
+
+    v1 = -v1;
+    sout | "-?:" | v1;
+
+    v1 = (vec4(float)){1.5f, 2.75f, -14.2f, -13.5f};
+    sout | "?+?:" | (v1 + (vec4(float)){0.8f, -0.3f, 5, 1});
+
+    v1 += (vec4(float)){0.8f, -0.3f, 5, 1};
+    sout | "?+=?:" | v1;
+
+    v1 = (vec4(float)){1.5f, 2.75f, 100.3f, -2.f};
+    sout | "v*s:" | v1 * 3.f;
+
+    sout | "s*v:" | 3.f * v1;
+
+    v1 *= 3;
+    sout | "?*=?:" | v1;
+
+    v1 = (vec4(float)){2, -0.1f, 45, -15};
+    sout | "?/?:" | (v1 / 3.f);
+
+    v1 /= 3.f;
+    sout | "?/=?:" | v1;
+
+    v1 = (vec4(float)){4,3,2,1};
+    v2 = (vec4(float)){0,4,-3,-6};
+    sout | "dot_1:" | dot(v1, v2);
+
+    v2 = (vec4(float)){1.3f, -2, 12.2f, 24.1345f};
+    sout | "dot_2:" | dot(v1, v2);
+
+    v1 = (vec4(float)){1,2,3,4};
+    sout | "length:" | length(v1);
+
+    sout | "length_squared:" | length_squared(v1);
+
+    v2 = (vec4(float)){6, -3.2f, 1, 12};
+    sout | "distance:" | distance(v1, v2);
+
+    sout | "normalize:" | normalize(v2);
+
+    sout | "project:" | project((vec4(float)){5,6,0,-1}, (vec4(float)){0, 0, 0, 1});
+    sout | "project_2:" | project((vec4(float)){5,6,3.2f,-16.75f}, (vec4(float)){10, 1, 4, 2});
+
+    v1 = (vec4(float)){5,6,2.333f,1};
+    v2 = (vec4(float)){1,0,-13.5f,2};
+    sout | "reflect:" | reflect(v1,v2);
+
+    v2 = (vec4(float)){0,-1,2,0.5f};
+    sout | "refract:" | refract(normalize(v1),normalize(v2),1.f);
+    sout | "refract:" | refract(normalize(v1),normalize(v2),1.f/1.33f);
+
+    vec4(float) geometric_normal = {5,6,1,2};
+    vec4(float) perturbed_normal = {4,5.5f,2,2};
+    vec4(float) eyeline = {-1,0.002f,-1.0345f,-2};
+    sout | "faceforward_nochange:" | faceforward(perturbed_normal, eyeline, geometric_normal);
+
+    eyeline = (vec4(float)){1,0.002f,-1.0345f,-2};
+    sout | "faceforward_flip:" | faceforward(perturbed_normal, eyeline, geometric_normal);
+}
Index: tests/vector_math/vec4_int.cfa
===================================================================
--- tests/vector_math/vec4_int.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
+++ tests/vector_math/vec4_int.cfa	(revision eef8dfb77fab582e377f7d45c1dca2dc31630bd9)
@@ -0,0 +1,130 @@
+#include <vec/vec4.hfa>
+#include <fstream.hfa>
+#include <limits.hfa>
+
+int main(void) {
+    vec4(int) v1 = {1,2,3,4};
+    sout | "ctor(x,y):" | v1;
+
+    vec4(int) v2 = v1;
+    sout | "copy ctor:" | v2;
+
+    {
+        vec4(int) u = {11, -5, 7, 11};
+        vec4(int) v = {2, 3, 7, 2};
+        sout | "?/?(vec)" | (u / v);
+        sout | "?/=?(vec)" | (u /= v);
+    }
+
+    {
+        vec4(int) u = {1, -3, 100, 1};
+        sout | "++?" | ++u;
+    }
+
+    {
+        vec4(int) u = {1, -3, 100, 3};
+        sout | "--?" | --u;
+    }
+    {
+        vec4(int) u = {1, -3, 100, 1};
+        sout | "?++(part 1)" | u++;
+        sout | "?++(part 2)" | u;
+    }
+    {
+        vec4(int) u = {1, -3, 100, 3};
+        sout | "?--(part 1)" | u--;
+        sout | "?--(part 2)" | u;
+    }
+
+    // %
+    {
+        vec4(int) u = {-12, 123, 100, 5};
+        sout | "?%?(int)" | (u % 2);
+        sout | "?%=?(int)" | (u %= 2);
+    }
+    {
+        vec4(int) u = {-12, 27, 100, 14};
+        vec4(int) v = {5, 4, 13, 5};
+        sout | "?%?(vec)" | (u % v);
+        sout | "?%=?(vec)" | (u %= v);
+    }
+
+
+    // &
+    {
+        vec4(int) u = {5, 1, 8, 15};
+        sout | "?&?(int)" | (u & (1 << 3));
+        sout | "?&=?(int)" | (u &= (1 << 3));
+    }
+    {
+        vec4(int) u = {5, 6, 5, 13};
+        vec4(int) v = {3, 7, 3, 7};
+        sout | "?&?(vec)" | (u & v);
+        sout | "?&=?(vec)" | (u &= v);
+    }
+
+    // |
+    {
+        vec4(int) u = {5, 1, 5, 8};
+        sout | "?|?(int)" | (u | 2);
+        sout | "?|=?(int)" | (u |= 2);
+    }
+    {
+        vec4(int) u = {5, 1, 5, 2};
+        vec4(int) v = {2, 1, 2, 8};
+        sout | "?|?(vec)" | (u | v);
+        sout | "?|=?(vec)" | (u |= v);
+    }
+
+    // ^
+    {
+        vec4(int) u = {9, 5, 8, 13};
+        sout | "?^?(int)" | (u ^ 5);
+        sout | "?^=?(int)" | (u ^= 5);
+    }
+    {
+        vec4(int) u = {9, 5, 8, 13};
+        vec4(int) v = {5, 0, 1, 7};
+        sout | "?^?(vec)" | (u ^ v);
+        sout | "?^=?(vec)" | (u ^= v);
+    }
+
+    // <<
+    {
+        vec4(int) u = {2, 1, 4, 6};
+        sout | "?<<?(int)" | (u << 1);
+        sout | "?<<=?(int)" | (u <<= 1);
+    }
+    {
+        vec4(int) u = {1, 1, 1, 4};
+        vec4(int) v = {5, 2, 1, 1};
+        sout | "?<<?(vec)" | (u << v);
+        sout | "?<<=?(vec)" | (u <<= v);
+    }
+
+    // >>
+    {
+        vec4(int) u = {2, 1, 4, 6};
+        sout | "?>>?(int)" | (u >> 1);
+        sout | "?>>=?(int)" | (u >>= 1);
+    }
+    {
+        vec4(int) u = {5, 7, 8, 32};
+        vec4(int) v = {1, 2, 2, 3};
+        sout | "?>>?(vec)" | (u >> v);
+        sout | "?>>=?(vec)" | (u >>= v);
+    }
+
+    // ~
+    {
+        vec4(int) u = {5, 2, 1, -1};
+        sout | "~?" | (~u);
+    }
+
+    v1 = (vec4(int)){2, 3, 4, 5};
+    v2 = (vec4(int)){13, 2, -7, -3};
+    sout | "dot:" | dot(v1, v2);
+
+    v1 = (vec4(int)){2,3,4,5};
+    sout | "length_squared:" | length_squared(v1);
+}
