Index: libcfa/src/Makefile.am
===================================================================
--- libcfa/src/Makefile.am	(revision 4f7b418b1407ea6b3c3af7a736618baceddcbbc1)
+++ libcfa/src/Makefile.am	(revision bdfc0321152dff26f371dfef5322533b5cdd57dc)
@@ -41,5 +41,6 @@
 headers_nosrc = math.hfa gmp.hfa time_t.hfa bits/align.hfa bits/containers.hfa bits/defs.hfa bits/debug.hfa bits/locks.hfa
 headers = fstream.hfa iostream.hfa iterator.hfa limits.hfa rational.hfa time.hfa stdlib.hfa common.hfa \
-	  containers/maybe.hfa containers/pair.hfa containers/result.hfa containers/vector.hfa
+	  containers/maybe.hfa containers/pair.hfa containers/result.hfa containers/vector.hfa \
+	  vec/vec.hfa vec/vec2.hfa vec/vec3.hfa vec/vec4.hfa
 
 libsrc = startup.cfa interpose.cfa bits/debug.cfa assert.cfa exception.c virtual.c heap.cfa ${headers:.hfa=.cfa}
Index: libcfa/src/Makefile.in
===================================================================
--- libcfa/src/Makefile.in	(revision 4f7b418b1407ea6b3c3af7a736618baceddcbbc1)
+++ libcfa/src/Makefile.in	(revision bdfc0321152dff26f371dfef5322533b5cdd57dc)
@@ -142,5 +142,6 @@
 	time.cfa stdlib.cfa common.cfa containers/maybe.cfa \
 	containers/pair.cfa containers/result.cfa \
-	containers/vector.cfa
+	containers/vector.cfa vec/vec.cfa vec/vec2.cfa vec/vec3.cfa \
+	vec/vec4.cfa
 am__dirstamp = $(am__leading_dot)dirstamp
 @BUILDLIB_TRUE@am__objects_1 = fstream.lo iostream.lo iterator.lo \
@@ -148,5 +149,6 @@
 @BUILDLIB_TRUE@	common.lo containers/maybe.lo \
 @BUILDLIB_TRUE@	containers/pair.lo containers/result.lo \
-@BUILDLIB_TRUE@	containers/vector.lo
+@BUILDLIB_TRUE@	containers/vector.lo vec/vec.lo vec/vec2.lo \
+@BUILDLIB_TRUE@	vec/vec3.lo vec/vec4.lo
 @BUILDLIB_TRUE@am__objects_2 = startup.lo interpose.lo bits/debug.lo \
 @BUILDLIB_TRUE@	assert.lo exception.lo virtual.lo heap.lo \
@@ -237,7 +239,8 @@
 	limits.hfa rational.hfa time.hfa stdlib.hfa common.hfa \
 	containers/maybe.hfa containers/pair.hfa containers/result.hfa \
-	containers/vector.hfa math.hfa gmp.hfa time_t.hfa \
-	bits/align.hfa bits/containers.hfa bits/defs.hfa \
-	bits/debug.hfa bits/locks.hfa concurrency/coroutine.hfa \
+	containers/vector.hfa vec/vec.hfa vec/vec2.hfa vec/vec3.hfa \
+	vec/vec4.hfa math.hfa gmp.hfa time_t.hfa bits/align.hfa \
+	bits/containers.hfa bits/defs.hfa bits/debug.hfa \
+	bits/locks.hfa concurrency/coroutine.hfa \
 	concurrency/thread.hfa concurrency/kernel.hfa \
 	concurrency/monitor.hfa concurrency/mutex.hfa \
@@ -460,5 +463,6 @@
 @BUILDLIB_FALSE@headers = 
 @BUILDLIB_TRUE@headers = fstream.hfa iostream.hfa iterator.hfa limits.hfa rational.hfa time.hfa stdlib.hfa common.hfa \
-@BUILDLIB_TRUE@	  containers/maybe.hfa containers/pair.hfa containers/result.hfa containers/vector.hfa
+@BUILDLIB_TRUE@	  containers/maybe.hfa containers/pair.hfa containers/result.hfa containers/vector.hfa \
+@BUILDLIB_TRUE@	  vec/vec.hfa vec/vec2.hfa vec/vec3.hfa vec/vec4.hfa
 
 @BUILDLIB_FALSE@libsrc = 
@@ -590,4 +594,14 @@
 containers/vector.lo: containers/$(am__dirstamp) \
 	containers/$(DEPDIR)/$(am__dirstamp)
+vec/$(am__dirstamp):
+	@$(MKDIR_P) vec
+	@: > vec/$(am__dirstamp)
+vec/$(DEPDIR)/$(am__dirstamp):
+	@$(MKDIR_P) vec/$(DEPDIR)
+	@: > vec/$(DEPDIR)/$(am__dirstamp)
+vec/vec.lo: vec/$(am__dirstamp) vec/$(DEPDIR)/$(am__dirstamp)
+vec/vec2.lo: vec/$(am__dirstamp) vec/$(DEPDIR)/$(am__dirstamp)
+vec/vec3.lo: vec/$(am__dirstamp) vec/$(DEPDIR)/$(am__dirstamp)
+vec/vec4.lo: vec/$(am__dirstamp) vec/$(DEPDIR)/$(am__dirstamp)
 
 libcfa.la: $(libcfa_la_OBJECTS) $(libcfa_la_DEPENDENCIES) $(EXTRA_libcfa_la_DEPENDENCIES) 
@@ -629,4 +643,6 @@
 	-rm -f containers/*.$(OBJEXT)
 	-rm -f containers/*.lo
+	-rm -f vec/*.$(OBJEXT)
+	-rm -f vec/*.lo
 
 distclean-compile:
@@ -694,4 +710,5 @@
 	-rm -rf concurrency/.libs concurrency/_libs
 	-rm -rf containers/.libs containers/_libs
+	-rm -rf vec/.libs vec/_libs
 install-nobase_cfa_includeHEADERS: $(nobase_cfa_include_HEADERS)
 	@$(NORMAL_INSTALL)
@@ -840,4 +857,6 @@
 	-rm -f containers/$(DEPDIR)/$(am__dirstamp)
 	-rm -f containers/$(am__dirstamp)
+	-rm -f vec/$(DEPDIR)/$(am__dirstamp)
+	-rm -f vec/$(am__dirstamp)
 
 maintainer-clean-generic:
Index: libcfa/src/vec/vec.hfa
===================================================================
--- libcfa/src/vec/vec.hfa	(revision bdfc0321152dff26f371dfef5322533b5cdd57dc)
+++ libcfa/src/vec/vec.hfa	(revision bdfc0321152dff26f371dfef5322533b5cdd57dc)
@@ -0,0 +1,120 @@
+#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 bdfc0321152dff26f371dfef5322533b5cdd57dc)
+++ libcfa/src/vec/vec2.hfa	(revision bdfc0321152dff26f371dfef5322533b5cdd57dc)
@@ -0,0 +1,273 @@
+#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 bdfc0321152dff26f371dfef5322533b5cdd57dc)
+++ libcfa/src/vec/vec3.hfa	(revision bdfc0321152dff26f371dfef5322533b5cdd57dc)
@@ -0,0 +1,282 @@
+#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 bdfc0321152dff26f371dfef5322533b5cdd57dc)
+++ libcfa/src/vec/vec4.hfa	(revision bdfc0321152dff26f371dfef5322533b5cdd57dc)
@@ -0,0 +1,278 @@
+#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);
+    }
+}
+