Index: libcfa/src/Makefile.am
===================================================================
--- libcfa/src/Makefile.am	(revision e112a2422e4e6bf9a9a4884df0f64f67835cc8ff)
+++ libcfa/src/Makefile.am	(revision d435c04cb8fd918cf352e0f7d2c50b5252396c26)
@@ -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 = \
@@ -94,4 +98,5 @@
 	concurrency/clib/cfathread.h \
 	concurrency/invoke.h \
+	concurrency/future.hfa \
 	concurrency/kernel/fwd.hfa
 
Index: libcfa/src/bits/locks.hfa
===================================================================
--- libcfa/src/bits/locks.hfa	(revision e112a2422e4e6bf9a9a4884df0f64f67835cc8ff)
+++ libcfa/src/bits/locks.hfa	(revision d435c04cb8fd918cf352e0f7d2c50b5252396c26)
@@ -283,4 +283,9 @@
 		void ^?{}(future_t &) {}
 
+		void reset(future_t & this) {
+			// needs to be in 0p or 1p
+			__atomic_exchange_n( &this.ptr, 0p, __ATOMIC_SEQ_CST);
+		}
+
 		// check if the future is available
 		bool available( future_t & this ) {
@@ -340,6 +345,12 @@
 
 		// Mark the future as abandoned, meaning it will be deleted by the server
-		void abandon( future_t & this ) {
+		bool abandon( future_t & this ) {
+			/* paranoid */ verify( this.ptr != 3p );
+
+			// Mark the future as abandonned
 			struct oneshot * got = __atomic_exchange_n( &this.ptr, 3p, __ATOMIC_SEQ_CST);
+
+			// If the future isn't already fulfilled, let the server delete it
+			if( got == 0p ) return false;
 
 			// got == 2p: the future is ready but the context hasn't fully been consumed
@@ -347,6 +358,11 @@
 			if( got == 2p ) {
 				while( this.ptr != 1p ) Pause();
-			}
-			return;
+				got = 1p;
+			}
+
+			// The future is completed delete it now
+			/* paranoid */ verify( this.ptr != 1p );
+			free( &this );
+			return true;
 		}
 
Index: libcfa/src/concurrency/future.hfa
===================================================================
--- libcfa/src/concurrency/future.hfa	(revision d435c04cb8fd918cf352e0f7d2c50b5252396c26)
+++ libcfa/src/concurrency/future.hfa	(revision d435c04cb8fd918cf352e0f7d2c50b5252396c26)
@@ -0,0 +1,119 @@
+//
+// Cforall Version 1.0.0 Copyright (C) 2020 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           : Thierry Delisle & Peiran Hong
+// Created On       : Wed Jan 06 17:33:18 2021
+// Last Modified By :
+// Last Modified On :
+// Update Count     :
+//
+
+#pragma once
+
+#include "bits/locks.hfa"
+#include "monitor.hfa"
+
+forall( otype T ) {
+	struct future {
+		inline future_t;
+		T result;
+	};
+
+	static inline {
+		// Reset future back to original state
+		void reset(future(T) & this) { reset( (future_t&)this ); }
+
+		// check if the future is available
+		bool available( future(T) & this ) { return available( (future_t&)this ); }
+
+		// Mark the future as abandoned, meaning it will be deleted by the server
+		// This doesn't work beause of the potential need for a destructor
+		void abandon( future(T) & this );
+
+		// Fulfil the future, returns whether or not someone was unblocked
+		bool fulfil( future(T) & this, T result ) {
+			this.result = result;
+			return fulfil( (future_t&)this );
+		}
+
+		// Wait for the future to be fulfilled
+		// Also return whether the thread had to block or not
+		[T, bool] wait( future(T) & this ) {
+			bool r = wait( (future_t&)this );
+			return [this.result, r];
+		}
+
+		// Wait for the future to be fulfilled
+		T wait( future(T) & this ) {
+			[T, bool] tt;
+			tt = wait(this);
+			return tt.0;
+		}
+	}
+}
+
+forall( otype T ) {
+	monitor multi_future {
+		inline future_t;
+		condition blocked;
+		bool has_first;
+		T result;
+	};
+
+	static inline {
+		void ?{}(multi_future(T) & this) {
+			this.has_first = false;
+		}
+
+		bool $first( multi_future(T) & mutex this ) {
+			if (this.has_first) {
+				wait( this.blocked );
+				return false;
+			}
+
+			this.has_first = true;
+			return true;
+		}
+
+		void $first_done( multi_future(T) & mutex this ) {
+			this.has_first = false;
+			signal_all( this.blocked );
+		}
+
+		// Reset future back to original state
+		void reset(multi_future(T) & mutex this) {
+			if( this.has_first != false) abort("Attempting to reset a multi_future with at least one blocked threads");
+			if( !is_empty(this.blocked) ) abort("Attempting to reset a multi_future with multiple blocked threads");
+			reset( (future_t&)this );
+		}
+
+		// Fulfil the future, returns whether or not someone was unblocked
+		bool fulfil( multi_future(T) & this, T result ) {
+			this.result = result;
+			return fulfil( (future_t&)this );
+		}
+
+		// Wait for the future to be fulfilled
+		// Also return whether the thread had to block or not
+		[T, bool] wait( multi_future(T) & this ) {
+			bool sw = $first( this );
+			bool w = !sw;
+			if ( sw ) {
+				w = wait( (future_t&)this );
+				$first_done( this );
+			}
+
+			return [this.result, w];
+		}
+
+		// Wait for the future to be fulfilled
+		T wait( multi_future(T) & this ) {
+			return wait(this).0;
+		}
+	}
+}
Index: libcfa/src/concurrency/monitor.hfa
===================================================================
--- libcfa/src/concurrency/monitor.hfa	(revision e112a2422e4e6bf9a9a4884df0f64f67835cc8ff)
+++ libcfa/src/concurrency/monitor.hfa	(revision d435c04cb8fd918cf352e0f7d2c50b5252396c26)
@@ -132,7 +132,8 @@
 
               void wait        ( condition & this, uintptr_t user_info = 0 );
+static inline bool is_empty    ( condition & this ) { return this.blocked.head == 1p; }
               bool signal      ( condition & this );
               bool signal_block( condition & this );
-static inline bool is_empty    ( condition & this ) { return this.blocked.head == 1p; }
+static inline bool signal_all  ( condition & this ) { bool ret = false; while(!is_empty(this)) { ret = signal(this) || ret; } return ret; }
          uintptr_t front       ( condition & this );
 
Index: libcfa/src/vec/vec.hfa
===================================================================
--- libcfa/src/vec/vec.hfa	(revision d435c04cb8fd918cf352e0f7d2c50b5252396c26)
+++ libcfa/src/vec/vec.hfa	(revision d435c04cb8fd918cf352e0f7d2c50b5252396c26)
@@ -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 d435c04cb8fd918cf352e0f7d2c50b5252396c26)
+++ libcfa/src/vec/vec2.hfa	(revision d435c04cb8fd918cf352e0f7d2c50b5252396c26)
@@ -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 d435c04cb8fd918cf352e0f7d2c50b5252396c26)
+++ libcfa/src/vec/vec3.hfa	(revision d435c04cb8fd918cf352e0f7d2c50b5252396c26)
@@ -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 d435c04cb8fd918cf352e0f7d2c50b5252396c26)
+++ libcfa/src/vec/vec4.hfa	(revision d435c04cb8fd918cf352e0f7d2c50b5252396c26)
@@ -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);
+    }
+}
+