diff --git a/.gitignore b/.gitignore
index 8862c08..e380892 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,4 @@
 *.[oa]
+*.exe
+test/libmatrix_test
+test/libmatrix_test.exe
diff --git a/Makefile b/Makefile
index 418311f..330c289 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,8 @@
-CXXFLAGS = -Wall -Werror -pedantic -O3
+COMMON_FLAGS = -std=gnu++26 -Wall -Werror -pedantic -march=native -O3
+ifeq ($(shell uname -m), x86_64)
+COMMON_FLAGS += -mfpmath=sse
+endif
+CXXFLAGS  ?= $(COMMON_FLAGS)
 LIBMATRIX = libmatrix.a
 LIBSRCS = mat.cc program.cc log.cc util.cc shader-source.cc
 LIBOBJS = $(LIBSRCS:.cc=.o)
diff --git a/log.cc b/log.cc
index 10603ae..9b5373b 100644
--- a/log.cc
+++ b/log.cc
@@ -10,7 +10,6 @@
 //     Alexandros Frantzis <alexandros.frantzis@linaro.org>
 //     Jesse Barker <jesse.barker@linaro.org>
 //
-#include <unistd.h>
 #include <cstdio>
 #include <cstdarg>
 #include <string>
@@ -22,6 +21,13 @@
 #include <android/log.h>
 #endif
 
+#ifdef _WIN32
+// On windows 'isatty' is found in <io.h>
+#include <io.h>
+#else
+#include <unistd.h>
+#endif
+
 using std::string;
 
 const string Log::continuation_prefix("\x10");
@@ -33,6 +39,7 @@ static const string terminal_color_normal("\033[0m");
 static const string terminal_color_red("\033[1;31m");
 static const string terminal_color_cyan("\033[36m");
 static const string terminal_color_yellow("\033[33m");
+static const string terminal_color_magenta("\033[35m");
 static const string empty;
 
 static void
@@ -96,6 +103,16 @@ print_prefixed_message(std::ostream& stream, const string& color, const string&
     delete[] buf;
 }
 
+#ifdef ANDROID
+static void
+android_vlog(int prio, const char *tag, const char *fmt, va_list ap)
+{
+    va_list aq;
+    va_copy(aq, ap);
+    __android_log_vprint(prio, tag, fmt, aq);
+    va_end(aq);
+}
+#endif
 
 void
 Log::info(const char *fmt, ...)
@@ -110,7 +127,7 @@ Log::info(const char *fmt, ...)
     const string& color(do_debug_ ? infocolor : empty);
     print_prefixed_message(std::cout, color, prefix, fmt, ap);
 #else
-    __android_log_vprint(ANDROID_LOG_INFO, appname_.c_str(), fmt, ap);
+    android_vlog(ANDROID_LOG_INFO, appname_.c_str(), fmt, ap);
 #endif
 
     if (extra_out_)
@@ -132,7 +149,7 @@ Log::debug(const char *fmt, ...)
     static const string& dbgcolor(isatty(fileno(stdout)) ? terminal_color_yellow : empty);
     print_prefixed_message(std::cout, dbgcolor, dbgprefix, fmt, ap);
 #else
-    __android_log_vprint(ANDROID_LOG_DEBUG, appname_.c_str(), fmt, ap);
+    android_vlog(ANDROID_LOG_DEBUG, appname_.c_str(), fmt, ap);
 #endif
 
     if (extra_out_)
@@ -152,7 +169,7 @@ Log::error(const char *fmt, ...)
     static const string& errcolor(isatty(fileno(stderr)) ? terminal_color_red : empty);
     print_prefixed_message(std::cerr, errcolor, errprefix, fmt, ap);
 #else
-    __android_log_vprint(ANDROID_LOG_ERROR, appname_.c_str(), fmt, ap);
+    android_vlog(ANDROID_LOG_ERROR, appname_.c_str(), fmt, ap);
 #endif
 
     if (extra_out_)
@@ -161,6 +178,26 @@ Log::error(const char *fmt, ...)
     va_end(ap);
 }
 
+void
+Log::warning(const char *fmt, ...)
+{
+    static const string warnprefix("Warning");
+    va_list ap;
+    va_start(ap, fmt);
+
+#ifndef ANDROID
+    static const string& warncolor(isatty(fileno(stderr)) ? terminal_color_magenta : empty);
+    print_prefixed_message(std::cerr, warncolor, warnprefix, fmt, ap);
+#else
+    android_vlog(ANDROID_LOG_WARN, appname_.c_str(), fmt, ap);
+#endif
+
+    if (extra_out_)
+        print_prefixed_message(*extra_out_, empty, warnprefix, fmt, ap);
+
+    va_end(ap);
+}
+
 void
 Log::flush()
 {
diff --git a/log.h b/log.h
index 9054323..573de78 100644
--- a/log.h
+++ b/log.h
@@ -32,6 +32,8 @@ class Log
     static void debug(const char *fmt, ...);
     // Emit an error message
     static void error(const char *fmt, ...);
+    // Emit a warning message
+    static void warning(const char *fmt, ...);
     // Explicit flush of the log buffer
     static void flush();
     // A prefix constant that informs the logging infrastructure that the log
diff --git a/mat.h b/mat.h
index 493af0c..f60bd17 100644
--- a/mat.h
+++ b/mat.h
@@ -15,6 +15,11 @@
 #include <iostream>
 #include <iomanip>
 #include "vec.h"
+#ifndef USE_EXCEPTIONS
+// If we're not throwing exceptions, we'll need the logger to make sure the
+// caller is informed of errors.
+#include "log.h"
+#endif // USE_EXCEPTIONS
 
 namespace LibMatrix
 {
@@ -103,12 +108,17 @@ class tmat2
     //
     // NOTE: If this is non-invertible, we will
     //       throw to avoid undefined behavior.
-    tmat2& inverse() throw(std::runtime_error)
+    tmat2& inverse()
     {
         T d(determinant());
         if (d == static_cast<T>(0))
         {
+#ifdef USE_EXCEPTIONS
             throw std::runtime_error("Matrix is noninvertible!!!!");
+#else // !USE_EXCEPTIONS
+            Log::error("Matrix is noninvertible!!!!\n");
+            return *this;
+#endif // USE_EXCEPTIONS
         }
         T c0r0(m_[3] / d);
         T c0r1(-m_[1] / d);
@@ -397,12 +407,17 @@ class tmat3
     //
     // NOTE: If this is non-invertible, we will
     //       throw to avoid undefined behavior.
-    tmat3& inverse() throw(std::runtime_error)
+    tmat3& inverse()
     {
         T d(determinant());
         if (d == static_cast<T>(0))
         {
+#ifdef USE_EXCEPTIONS
             throw std::runtime_error("Matrix is noninvertible!!!!");
+#else // !USE_EXCEPTIONS
+            Log::error("Matrix is noninvertible!!!!\n");
+            return *this;
+#endif // USE_EXCEPTIONS
         }
         tmat2<T> minor0(m_[4], m_[5], m_[7], m_[8]);
         tmat2<T> minor1(m_[7], m_[8], m_[1], m_[2]);
@@ -771,12 +786,17 @@ class tmat4
     //
     // NOTE: If this is non-invertible, we will
     //       throw to avoid undefined behavior.
-    tmat4& inverse() throw(std::runtime_error)
+    tmat4& inverse()
     {
         T d(determinant());
         if (d == static_cast<T>(0))
         {
+#ifdef USE_EXCEPTIONS
             throw std::runtime_error("Matrix is noninvertible!!!!");
+#else // !USE_EXCEPTIONS
+            Log::error("Matrix is noninvertible!!!!\n");
+            return *this;
+#endif // USE_EXCEPTIONS
         }
         tmat3<T> minor0(m_[5], m_[6], m_[7], m_[9], m_[10], m_[11], m_[13], m_[14], m_[15]);
         tmat3<T> minor1(m_[1], m_[2], m_[3], m_[13], m_[14], m_[15], m_[9], m_[10], m_[11]);
diff --git a/program.cc b/program.cc
index b27298b..c952d60 100644
--- a/program.cc
+++ b/program.cc
@@ -190,7 +190,7 @@ Program::addShader(unsigned int type, const string& source)
     }
 
     shader.attach(handle_);
-    shaders_.push_back(shader);
+    shaders_.push_back(std::move(shader));
     return;
 }
 
diff --git a/program.h b/program.h
index b7ba3df..1c9de08 100644
--- a/program.h
+++ b/program.h
@@ -15,6 +15,7 @@
 #include <string>
 #include <vector>
 #include <map>
+#include <utility>
 #include "mat.h"
 
 // Simple shader container.  Abstracts all of the OpenGL bits, but leaves
@@ -35,6 +36,13 @@ class Shader
         message_(shader.message_),
         ready_(shader.ready_),
         valid_(shader.valid_) {}
+    Shader(Shader&& shader) noexcept:
+        handle_(std::exchange(shader.handle_, 0)),
+        type_(std::exchange(shader.type_, 0)),
+        source_(std::move(shader.source_)),
+        message_(std::move(shader.message_)),
+        ready_(std::exchange(shader.ready_, false)),
+        valid_(std::exchange(shader.valid_, false)) {}
     Shader(unsigned int type, const std::string& source);
     ~Shader();
 
diff --git a/shader-source.cc b/shader-source.cc
index d34c8a0..659832e 100644
--- a/shader-source.cc
+++ b/shader-source.cc
@@ -18,6 +18,36 @@
 #include "vec.h"
 #include "util.h"
 
+namespace
+{
+
+bool is_valid_precision_value(ShaderSource::PrecisionValue precision_value)
+{
+    switch(precision_value) {
+        case ShaderSource::PrecisionValueLow:
+        case ShaderSource::PrecisionValueMedium:
+        case ShaderSource::PrecisionValueHigh:
+        case ShaderSource::PrecisionValueDefault:
+            return true;
+        default:
+            return false;
+    }
+}
+
+bool is_valid_shader_type(ShaderSource::ShaderType shader_type)
+{
+    switch(shader_type) {
+        case ShaderSource::ShaderTypeVertex:
+        case ShaderSource::ShaderTypeFragment:
+        case ShaderSource::ShaderTypeUnknown:
+            return true;
+        default:
+            return false;
+    }
+}
+
+}
+
 /**
  * Holds default precision values for all shader types
  * (even the unknown type, which is hardwired to default precision values)
@@ -34,7 +64,7 @@ ShaderSource::default_precision_(ShaderSource::ShaderTypeUnknown + 1);
 bool
 ShaderSource::load_file(const std::string& filename, std::string& str)
 {
-    std::auto_ptr<std::istream> is_ptr(Util::get_resource(filename));
+    std::unique_ptr<std::istream> is_ptr(Util::get_resource(filename));
     std::istream& inputFile(*is_ptr);
 
     if (!inputFile)
@@ -421,7 +451,9 @@ ShaderSource::emit_precision(std::stringstream& ss, ShaderSource::PrecisionValue
             ss << "#endif" << std::endl;
         }
     }
-    else if (val >= 0 && val < ShaderSource::PrecisionValueDefault) {
+    else if (is_valid_precision_value(val) &&
+             val != ShaderSource::PrecisionValueDefault)
+    {
         ss << "precision " << precision_map[val] << " ";
         ss << type_str << ";" << std::endl;
     }
@@ -456,6 +488,22 @@ ShaderSource::str()
         precision = default_precision(type_);
 
     /* Create the precision statements */
+    std::stringstream precision_macros_ss;
+
+    precision_macros_ss << "#if defined(GL_ES)";
+    if (type_ == ShaderSource::ShaderTypeFragment)
+        precision_macros_ss << " && defined(GL_FRAGMENT_PRECISION_HIGH)";
+    precision_macros_ss << std::endl;
+    precision_macros_ss << "#define HIGHP_OR_DEFAULT highp" << std::endl;
+    precision_macros_ss << "#else" << std::endl;
+    precision_macros_ss << "#define HIGHP_OR_DEFAULT" << std::endl;
+    precision_macros_ss << "#endif" << std::endl;
+    precision_macros_ss << "#if defined(GL_ES)"  << std::endl;
+    precision_macros_ss << "#define MEDIUMP_OR_DEFAULT mediump" << std::endl;
+    precision_macros_ss << "#else" << std::endl;
+    precision_macros_ss << "#define MEDIUMP_OR_DEFAULT" << std::endl;
+    precision_macros_ss << "#endif" << std::endl;
+
     std::stringstream ss;
 
     emit_precision(ss, precision.int_precision, "int");
@@ -469,7 +517,7 @@ ShaderSource::str()
         precision_str.insert(precision_str.size(), "#endif\n");
     }
 
-    return precision_str + source_.str();
+    return precision_macros_ss.str() + precision_str + source_.str();
 }
 
 /**
@@ -512,7 +560,7 @@ void
 ShaderSource::default_precision(const ShaderSource::Precision& precision,
                                 ShaderSource::ShaderType type)
 {
-    if (type < 0 || type > ShaderSource::ShaderTypeUnknown)
+    if (!is_valid_shader_type(type))
         type = ShaderSource::ShaderTypeUnknown;
 
     if (type == ShaderSource::ShaderTypeUnknown) {
@@ -537,7 +585,7 @@ ShaderSource::default_precision(const ShaderSource::Precision& precision,
 const ShaderSource::Precision&
 ShaderSource::default_precision(ShaderSource::ShaderType type)
 {
-    if (type < 0 || type > ShaderSource::ShaderTypeUnknown)
+    if (!is_valid_shader_type(type))
         type = ShaderSource::ShaderTypeUnknown;
 
     return default_precision_[type];
diff --git a/util.cc b/util.cc
index d96f393..c113a86 100644
--- a/util.cc
+++ b/util.cc
@@ -12,11 +12,15 @@
 //
 #include <sstream>
 #include <fstream>
-#include <sys/time.h>
+#include <chrono>
 #ifdef ANDROID
 #include <android/asset_manager.h>
+#endif
+#ifdef _WIN32
+#include <windows.h>
 #else
-#include <dirent.h>
+#include <unistd.h>
+#include <sys/resource.h>
 #endif
 
 #include "log.h"
@@ -101,6 +105,7 @@ fill_escape_vector(const string &str, vector<bool> &esc_vec)
                     state = StateNormal;
                 else
                     esc = true;
+                break;
             default:
                 break;
         }
@@ -220,58 +225,33 @@ Util::split(const string& src, char delim, vector<string>& elementVec,
 uint64_t
 Util::get_timestamp_us()
 {
-    struct timeval tv;
-    gettimeofday(&tv, NULL);
-    uint64_t now = static_cast<uint64_t>(tv.tv_sec) * 1000000 +
-                   static_cast<double>(tv.tv_usec);
-    return now;
-}
-
-std::string
-Util::appname_from_path(const std::string& path)
-{
-    std::string::size_type slashPos = path.rfind("/");
-    std::string::size_type startPos(0);
-    if (slashPos != std::string::npos)
-    {
-        startPos = slashPos + 1;
-    }
-    return std::string(path, startPos, std::string::npos);
+    return
+        std::chrono::duration_cast<std::chrono::microseconds>(
+            std::chrono::steady_clock::now().time_since_epoch()).count();
 }
 
 #ifndef ANDROID
 
 std::istream *
-Util::get_resource(const std::string &path)
+Util::get_resource(const std::filesystem::path &path)
 {
-    std::ifstream *ifs = new std::ifstream(path.c_str());
+    std::ifstream *ifs = new std::ifstream(path, std::ios::binary);
 
     return static_cast<std::istream *>(ifs);
 }
 
 void
-Util::list_files(const std::string& dirName, std::vector<std::string>& fileVec)
+Util::list_files(const std::filesystem::path& dirName,
+                 std::vector<std::filesystem::path>& fileVec)
 {
-    DIR* dir = opendir(dirName.c_str());
-    if (!dir)
+    try
     {
-        Log::error("Failed to open models directory '%s'\n", dirName.c_str());
-        return;
+        for (const auto& entry : std::filesystem::directory_iterator{dirName})
+            fileVec.push_back(entry.path());
     }
-
-    struct dirent* entry = readdir(dir);
-    while (entry)
+    catch (...)
     {
-        std::string pathname(dirName + "/");
-        pathname += std::string(entry->d_name);
-        // Skip '.' and '..'
-        if (entry->d_name[0] != '.')
-        {
-            fileVec.push_back(pathname);
-        }
-        entry = readdir(dir);
     }
-    closedir(dir);
 }
 
 #else
@@ -291,9 +271,9 @@ Util::android_get_asset_manager()
 }
 
 std::istream *
-Util::get_resource(const std::string &path)
+Util::get_resource(const std::filesystem::path &path)
 {
-    std::string path2(path);
+    std::string path2 = path.string();
     /* Remove leading '/' from path name, it confuses the AssetManager */
     if (path2.size() > 0 && path2[0] == '/')
         path2.erase(0, 1);
@@ -315,10 +295,11 @@ Util::get_resource(const std::string &path)
 }
 
 void
-Util::list_files(const std::string& dirName, std::vector<std::string>& fileVec)
+Util::list_files(const std::filesystem::path& dirName,
+                 std::vector<std::filesystem::path>& fileVec)
 {
     AAssetManager *mgr(Util::android_get_asset_manager());
-    std::string dir_name(dirName);
+    std::string dir_name = dirName.string();
 
     /* Remove leading '/' from path, it confuses the AssetManager */
     if (dir_name.size() > 0 && dir_name[0] == '/')
@@ -341,3 +322,60 @@ Util::list_files(const std::string& dirName, std::vector<std::string>& fileVec)
     AAssetDir_close(dir);
 }
 #endif
+
+unsigned int
+Util::get_num_processors()
+{
+#ifdef _WIN32
+    SYSTEM_INFO sysinfo;
+    GetSystemInfo(&sysinfo);
+    return sysinfo.dwNumberOfProcessors;
+#else
+    return sysconf(_SC_NPROCESSORS_ONLN);
+#endif
+}
+
+void
+Util::get_process_times(double *user_sec, double *system_sec)
+{
+#ifdef _WIN32
+    FILETIME creationTime, exitTime, kernelTime, userTime;
+    ULARGE_INTEGER user, kernel;
+
+    GetProcessTimes(GetCurrentProcess(), &creationTime, &exitTime, &kernelTime, &userTime);
+    // convert FILETIME to ULARGE_INTEGER
+    user.LowPart = userTime.dwLowDateTime;
+    user.HighPart = userTime.dwHighDateTime;
+    kernel.LowPart = kernelTime.dwLowDateTime;
+    kernel.HighPart = kernelTime.dwHighDateTime;
+    // FILETIME contains the number of 100 nsec intervals.
+    *user_sec = user.QuadPart / 1e7;
+    *system_sec = kernel.QuadPart / 1e7;
+#else
+    struct rusage usage;
+
+    getrusage(RUSAGE_SELF, &usage);
+    *user_sec = usage.ru_utime.tv_sec + usage.ru_utime.tv_usec / 1e6;
+    *system_sec = usage.ru_stime.tv_sec + usage.ru_stime.tv_usec / 1e6;
+#endif
+}
+
+double
+Util::get_idle_time()
+{
+#ifdef _WIN32
+    FILETIME idleTime;
+    GetSystemTimes(&idleTime, nullptr, nullptr);
+    ULARGE_INTEGER ulIdleTime;
+    ulIdleTime.LowPart = idleTime.dwLowDateTime;
+    ulIdleTime.HighPart = idleTime.dwHighDateTime;
+    // FILETIME contains the number of 100 nsec intervals.
+    return ulIdleTime.QuadPart / 1e7;
+#else
+    double uptime, idle;
+    std::ifstream ifs("/proc/uptime");
+    ifs >> uptime >> idle;
+    if (!ifs.fail()) return idle;
+    return 0.0;
+#endif
+}
diff --git a/util.h b/util.h
index 2b0f0f0..bef592e 100644
--- a/util.h
+++ b/util.h
@@ -16,7 +16,9 @@
 #include <string>
 #include <vector>
 #include <istream>
+#include <iomanip>
 #include <sstream>
+#include <filesystem>
 #include <stdint.h>
 
 #ifdef ANDROID
@@ -64,7 +66,7 @@ struct Util {
      * Returns a pointer to an input stream, which must be deleted when no
      * longer in use.
      */
-    static std::istream *get_resource(const std::string &path);
+    static std::istream *get_resource(const std::filesystem::path &path);
     /**
      * list_files() - Get a list of the files in a given directory.
      *
@@ -74,7 +76,8 @@ struct Util {
      * Obtains a list of the files in @dirName, and returns them in the string
      * vector @fileVec.
      */
-    static void list_files(const std::string& dirName, std::vector<std::string>& fileVec);
+    static void list_files(const std::filesystem::path& dirName,
+                           std::vector<std::filesystem::path>& fileVec);
     /**
      * dispose_pointer_vector() - cleans up a vector of pointers
      *
@@ -105,7 +108,7 @@ struct Util {
     {
         std::stringstream ss(asString);
         T retVal = T();
-        ss >> retVal;
+        ss >> std::setbase(0) >> retVal;
         return retVal;
     }
     /**
@@ -122,14 +125,22 @@ struct Util {
         return ss.str();
     }
     /**
-     * appname_from_path() - get the name of an executable from an absolute path
+     * toString() - Converts a double type to a string with precision.
      *
-     * @path:   absolute path of the running application (argv[0])
-     *
-     * Returns the last portion of @path (everything after the final '/').
+     * @t:         a double value to be converted to a string
+     * @precision: the precision to use for the conversion
      */
     static std::string
-    appname_from_path(const std::string& path);
+    toString(double t, int precision)
+    {
+        std::stringstream ss;
+        ss << std::fixed << std::setprecision(precision) << t;
+        return ss.str();
+    }
+
+    static unsigned int get_num_processors();
+    static void get_process_times(double *user_sec, double *system_sec);
+    static double get_idle_time();
 
 #ifdef ANDROID
     static void android_set_asset_manager(AAssetManager *asset_manager);
diff --git a/vec.h b/vec.h
index 2680ebc..1797087 100644
--- a/vec.h
+++ b/vec.h
@@ -14,392 +14,157 @@
 
 #include <iostream> // only needed for print() functions...
 #include <math.h>
+#include <array>
+#include <bit>
+#include <numeric>
+#include <utility>
 
-namespace LibMatrix
-{
-// A template class for creating, managing and operating on a 2-element vector
-// of any type you like (intended for built-in types, but as long as it 
-// supports the basic arithmetic and assignment operators, any type should
-// work).
 template<typename T>
-class tvec2
-{
-public:
-    tvec2() :
-        x_(0),
-        y_(0) {}
-    tvec2(const T t) :
-        x_(t),
-        y_(t) {}
-    tvec2(const T x, const T y) :
-        x_(x),
-        y_(y) {}
-    tvec2(const tvec2& v) :
-        x_(v.x_),
-        y_(v.y_) {}
-    ~tvec2() {}
-
-    // Print the elements of the vector to standard out.
-    // Really only useful for debug and test.
-    void print() const
-    {
-        std::cout << "| " << x_ << " " << y_ << " |" << std::endl;
-    }
-
-    // Allow raw data access for API calls and the like.
-    // For example, it is valid to pass a tvec2<float> into a call to
-    // the OpenGL command "glUniform2fv()".
-    operator const T*() const { return &x_;}
-
-    // Get and set access members for the individual elements.
-    const T x() const { return x_; }
-    const T y() const { return y_; }
-
-    void x(const T& val) { x_ = val; }
-    void y(const T& val) { y_ = val; }
-
-    // A direct assignment of 'rhs' to this.  Return a reference to this.
-    tvec2& operator=(const tvec2& rhs)
-    {
-        if (this != &rhs)
-        {
-            x_ = rhs.x_;
-            y_ = rhs.y_;
-        }
-        return *this;
-    }
-
-    // Divide this by a scalar.  Return a reference to this.
-    tvec2& operator/=(const T& rhs)
-    {
-        x_ /= rhs;
-        y_ /= rhs;
-        return *this;
-    }
-
-    // Divide a copy of this by a scalar.  Return the copy.
-    const tvec2 operator/(const T& rhs) const
-    {
-        return tvec2(*this) /= rhs;
-    }
-
-    // Component-wise divide of this by another vector.
-    // Return a reference to this.
-    tvec2& operator/=(const tvec2& rhs)
-    {
-        x_ /= rhs.x_;
-        y_ /= rhs.y_;
-        return *this;
-    }
-
-    // Component-wise divide of a copy of this by another vector.
-    // Return the copy.
-    const tvec2 operator/(const tvec2& rhs) const
-    {
-        return tvec2(*this) /= rhs;
-    }
-
-    // Multiply this by a scalar.  Return a reference to this.
-    tvec2& operator*=(const T& rhs)
-    {
-        x_ *= rhs;
-        y_ *= rhs;
-        return *this;
-    }
-
-    // Multiply a copy of this by a scalar.  Return the copy.
-    const tvec2 operator*(const T& rhs) const
-    {
-        return tvec2(*this) *= rhs;
-    }
-
-    // Component-wise multiply of this by another vector.
-    // Return a reference to this.
-    tvec2& operator*=(const tvec2& rhs)
-    {
-        x_ *= rhs.x_;
-        y_ *= rhs.y_;
-        return *this;
-    }
-
-    // Component-wise multiply of a copy of this by another vector.
-    // Return the copy.
-    const tvec2 operator*(const tvec2& rhs) const
-    {
-        return tvec2(*this) *= rhs;
-    }
-
-    // Add a scalar to this.  Return a reference to this.
-    tvec2& operator+=(const T& rhs)
-    {
-        x_ += rhs;
-        y_ += rhs;
-        return *this;
-    }
-    
-    // Add a scalar to a copy of this.  Return the copy.
-    const tvec2 operator+(const T& rhs) const
-    {
-        return tvec2(*this) += rhs;
-    }
-
-    // Component-wise addition of another vector to this.
-    // Return a reference to this.
-    tvec2& operator+=(const tvec2& rhs)
-    {
-        x_ += rhs.x_;
-        y_ += rhs.y_;
-        return *this;
-    }
-
-    // Component-wise addition of another vector to a copy of this.
-    // Return the copy.
-    const tvec2 operator+(const tvec2& rhs) const
-    {
-        return tvec2(*this) += rhs;
-    }
-
-    // Subtract a scalar from this.  Return a reference to this.
-    tvec2& operator-=(const T& rhs)
-    {
-        x_ -= rhs;
-        y_ -= rhs;
-        return *this;
-    }
-    
-    // Subtract a scalar from a copy of this.  Return the copy.
-    const tvec2 operator-(const T& rhs) const
-    {
-        return tvec2(*this) -= rhs;
-    }
-
-    // Component-wise subtraction of another vector from this.
-    // Return a reference to this.
-    tvec2& operator-=(const tvec2& rhs)
-    {
-        x_ -= rhs.x_;
-        y_ -= rhs.y_;
-        return *this;
-    }
+concept scalar = std::integral<T> || std::floating_point<T>;
 
-    // Component-wise subtraction of another vector from a copy of this.
-    // Return the copy.
-    const tvec2 operator-(const tvec2& rhs) const
-    {
-        return tvec2(*this) -= rhs;
-    }
-
-    // Compute the length of this and return it.
-    float length() const
-    {
-        return sqrt(dot(*this, *this));
-    }
+template<typename T>
+concept fscalar = std::floating_point<T>;
 
-    // Make this a unit vector.
-    void normalize()
-    {
-        float l = length();
-        x_ /= l;
-        y_ /= l;
-    }
+template<typename T>
+concept iscalar = std::integral<T>;
 
-    // Compute the dot product of two vectors.
-    static T dot(const tvec2& v1, const tvec2& v2)
-    {
-        return (v1.x_ * v2.x_) + (v1.y_ * v2.y_); 
-    }
+template<typename T>
+concept uscalar = std::unsigned_integral<T>;
 
-private:
-    T x_;
-    T y_;
+namespace LibMatrix
+{
+enum class align : size_t
+{
+	none     = 0,
+	element  = 1,
+	vector   = 2,
+	adaptive = 3
 };
 
-// A template class for creating, managing and operating on a 3-element vector
-// of any type you like (intended for built-in types, but as long as it 
-// supports the basic arithmetic and assignment operators, any type should
-// work).
-template<typename T>
-class tvec3
+// aligned n-element vector based on std:array compatible with every kind of SIMD optimization
+template<scalar T, size_t N = 3, enum align A = align::adaptive, size_t N_POW2 = std::bit_ceil<size_t>(N)>
+struct alignas(((N == N_POW2 && A != align::element) || A == align::vector) ? N_POW2 * sizeof(T) : sizeof(T)) tvec : std::array<T,N>
 {
-public:
-    tvec3() :
-        x_(0),
-        y_(0),
-        z_(0) {}
-    tvec3(const T t) :
-        x_(t),
-        y_(t),
-        z_(t) {}
-    tvec3(const T x, const T y, const T z) :
-        x_(x),
-        y_(y),
-        z_(z) {}
-    tvec3(const tvec3& v) :
-        x_(v.x_),
-        y_(v.y_),
-        z_(v.z_) {}
-    ~tvec3() {}
-
-    // Print the elements of the vector to standard out.
-    // Really only useful for debug and test.
-    void print() const
-    {
-        std::cout << "| " << x_ << " " << y_ << " " << z_ << " |" << std::endl;
-    }
+    tvec() { (*this).fill((T)0); }
+    tvec(const T& t) { (*this).fill((T)t); }
 
-    // Allow raw data access for API calls and the like.
-    // For example, it is valid to pass a tvec3<float> into a call to
-    // the OpenGL command "glUniform3fv()".
-    operator const T*() const { return &x_;}
+    template<scalar... I> requires((sizeof...(I) > 1) && (sizeof...(I) <= N))
+    tvec(const I... args) : std::array<T,N>{{ (T)args... }} {}
 
-    // Get and set access members for the individual elements.
-    const T x() const { return x_; }
-    const T y() const { return y_; }
-    const T z() const { return z_; }
+    template<size_t N_LHS, enum align A_RHS = align::adaptive> 
+    operator tvec<T,N_LHS,A_RHS>() { return *reinterpret_cast<tvec<T,N_LHS,A_RHS>*>(this); }
+    template<size_t N_LHS, enum align A_RHS = align::adaptive> 
+    operator const tvec<T,N_LHS,A_RHS>() const { return *reinterpret_cast<const tvec<T,N_LHS,A_RHS>*>(this); }
 
-    void x(const T& val) { x_ = val; }
-    void y(const T& val) { y_ = val; }
-    void z(const T& val) { z_ = val; }
+    template<enum align A_RHS = align::adaptive>
+    tvec(const tvec<T,3>& src, const T w = 1) requires (N > 3) { (*this).fill((T)0); (*this)[0] = src[0]; (*this)[1] = src[1]; (*this)[2] = src[2]; (*this)[3] = w; };
 
-    // A direct assignment of 'rhs' to this.  Return a reference to this.
-    tvec3& operator=(const tvec3& rhs)
+    void print() const
     {
-        if (this != &rhs)
-        {
-            x_ = rhs.x_;
-            y_ = rhs.y_;
-            z_ = rhs.z_;
-        }
-        return *this;
+        std::cout << "| ";
+        for(T& i : (*this))
+		std::cout << i << " ";
+	std::cout << "|" << std::endl;
     }
 
-    // Divide this by a scalar.  Return a reference to this.
-    tvec3& operator/=(const T& rhs)
-    {
-        x_ /= rhs;
-        y_ /= rhs;
-        z_ /= rhs;
-        return *this;
-    }
+    operator const T*() const { return (*this).data(); }
 
-    // Divide a copy of this by a scalar.  Return the copy.
-    const tvec3 operator/(const T& rhs) const
-    {
-        return tvec3(*this) /= rhs;
-    }
+    // Get and set access members for the individual elements.
+    const T x() const                 { return (*this)[0]; }
+    const T y() const requires(N > 1) { return (*this)[1]; }
+    const T z() const requires(N > 2) { return (*this)[2]; }
+    const T w() const requires(N > 3) { return (*this)[3]; }
 
-    // Component-wise divide of this by another vector.
-    // Return a reference to this.
-    tvec3& operator/=(const tvec3& rhs)
-    {
-        x_ /= rhs.x_;
-        y_ /= rhs.y_;
-        z_ /= rhs.z_;
-        return *this;
-    }
+    void x(const T& val)                 { (*this)[0] = val; }
+    void y(const T& val) requires(N > 1) { (*this)[1] = val; }
+    void z(const T& val) requires(N > 2) { (*this)[2] = val; }
+    void w(const T& val) requires(N > 3) { (*this)[3] = val; }
 
-    // Component-wise divide of a copy of this by another vector.
-    // Return the copy.
-    const tvec3 operator/(const tvec3& rhs) const
-    {
-        return tvec3(*this) /= rhs;
-    }
+    const tvec<T,3,A> yzx() const requires(N > 2) { return { (*this)[1], (*this)[2], (*this)[0] }; }
+    const tvec<T,3,A> zxy() const requires(N > 2) { return reinterpret_cast<tvec<T,3,A> const &>(*this); }
+    const tvec<T,3,A> xyz(T w = 1) requires (N == 4) { w *= w(); return (w != (T)0 && w != (T)1) ? reinterpret_cast<tvec<T,3,A> const &>(*this) / w : reinterpret_cast<tvec<T,3,A> const &>(*this); }
+    const tvec<T,3,A> yxz() const requires(N > 2) { return { (*this)[1], (*this)[0], (*this)[2] }; }
+    const tvec<T,4,A> xyzw() const requires(N > 3) { return reinterpret_cast<tvec<T,4,A> const &>(*this); }
 
-    // Multiply this by a scalar.  Return a reference to this.
-    tvec3& operator*=(const T& rhs)
-    {
-        x_ *= rhs;
-        y_ *= rhs;
-        z_ *= rhs;
-        return *this;
-    }
+    const tvec<T,4,A> position() requires(N > 2) { return tvec<T,4,A>(xyz(), 1); }
+    const tvec<T,4,A> direction() requires(N > 2) { return tvec<T,4,A>(xyz(), 0); }
 
-    // Multiply a copy of this by a scalar.  Return the copy.
-    const tvec3 operator*(const T& rhs) const
+    template<scalar T_RHS = T, size_t N_RHS = N, enum align A_RHS = A, scalar T_DST = decltype((T)1 * (T_RHS)1)>
+    inline constexpr const tvec<T,N,A> operator*(const tvec<T_RHS,N_RHS,A_RHS>& rhs) const
     {
-        return tvec3(*this) *= rhs;
+        tvec<T, std::min((*this).size(), rhs.size())> dst = {};
+        std::transform((*this).cbegin(),(*this).cbegin() + dst.size(), rhs.cbegin(), dst.begin(), std::multiplies<>{});
+	return dst;
     }
-
-    // Component-wise multiply of this by another vector.
-    // Return a reference to this.
-    tvec3& operator*=(const tvec3& rhs)
+    template<scalar T_RHS = T, size_t N_RHS = N, enum align A_RHS = A, scalar T_DST = decltype((T)1 / (T_RHS)1)>
+    inline constexpr const tvec<T,N,A> operator/(const tvec<T_RHS,N_RHS,A_RHS>& rhs) const
     {
-        x_ *= rhs.x_;
-        y_ *= rhs.y_;
-        z_ *= rhs.z_;
-        return *this;
+        tvec<T, std::min((*this).size(), rhs.size())> dst = {};
+	std::transform((*this).cbegin(),(*this).cbegin() + dst.size(), rhs.cbegin(), dst.begin(), std::divides<>{});
+	return dst;
     }
-
-    // Component-wise multiply of a copy of this by another vector.
-    // Return the copy.
-    const tvec3 operator*(const tvec3& rhs) const
+    template<scalar T_RHS = T, size_t N_RHS = N, enum align A_RHS = A, scalar T_DST = decltype((T)1 + (T_RHS)1)>
+    inline constexpr const tvec<T_DST,N,A> operator+(const tvec<T_RHS,N_RHS,A_RHS>& rhs) const
     {
-        return tvec3(*this) *= rhs;
+        tvec<T, std::min((*this).size(), rhs.size())> dst = {};
+	std::transform((*this).cbegin(),(*this).cbegin() + dst.size(), rhs.cbegin(), dst.begin(), std::plus<>{});
+	return dst;
     }
-
-    // Add a scalar to this.  Return a reference to this.
-    tvec3& operator+=(const T& rhs)
+    template<scalar T_RHS = T, size_t N_RHS = N, enum align A_RHS = A, scalar T_DST = decltype((T)1 - (T_RHS)1)>
+    inline constexpr const tvec<T_DST,N,A> operator-(const tvec<T_RHS,N_RHS,A_RHS>& rhs) const
     {
-        x_ += rhs;
-        y_ += rhs;
-        z_ += rhs;
-        return *this;
+        tvec<T, std::min((*this).size(), rhs.size())> dst = {};
+	std::transform((*this).cbegin(),(*this).cbegin() + dst.size(), rhs.cbegin(), dst.begin(), std::minus<>{});
+	return dst;
     }
-
-    // Add a scalar to a copy of this.  Return the copy.
-    const tvec3 operator+(const T& rhs) const
+    /* arithmetic scalar operators with constructor fill */
+    template<scalar T_RHS = T>
+    inline constexpr const tvec<T,N,A> operator*(const T_RHS& rhs) const
     {
-        return tvec3(*this) += rhs;
+        tvec<T, N, A> dst((T)rhs);
+	std::transform((*this).cbegin(),(*this).cend(), dst.cbegin(), dst.begin(), std::multiplies<>{});
+	return dst;
     }
-
-    // Component-wise addition of another vector to this.
-    // Return a reference to this.
-    tvec3& operator+=(const tvec3& rhs)
+    template<scalar T_RHS = T>
+    inline constexpr const tvec<T,N,A> operator/(const T_RHS& rhs) const
     {
-        x_ += rhs.x_;
-        y_ += rhs.y_;
-        z_ += rhs.z_;
-        return *this;
+        tvec<T, N, A> dst((T)rhs);
+	std::transform((*this).cbegin(),(*this).cend(), dst.cbegin(), dst.begin(), std::divides<>{});
+	return dst;
     }
-
-    // Component-wise addition of another vector to a copy of this.
-    // Return the copy.
-    const tvec3 operator+(const tvec3& rhs) const
+    template<scalar T_RHS = T>
+    inline constexpr const tvec<T,N,A> operator+(const T_RHS& rhs) const
     {
-        return tvec3(*this) += rhs;
+        tvec<T, N, A> dst((T)rhs);
+	std::transform((*this).cbegin(),(*this).cend(), dst.cbegin(), dst.begin(), std::plus<>{});
+	return dst;
     }
-
-    // Subtract a scalar from this.  Return a reference to this.
-    tvec3& operator-=(const T& rhs)
+    template<scalar T_RHS = T>
+    inline constexpr const tvec<T,N,A> operator-(const T_RHS& rhs) const
     {
-        x_ -= rhs;
-        y_ -= rhs;
-        z_ -= rhs;
-        return *this;
+        tvec<T, N, A> dst((T)rhs);
+	std::transform((*this).cbegin(),(*this).cend(), dst.cbegin(), dst.begin(), std::minus<>{});
+	return dst;
     }
 
-    // Subtract a scalar from a copy of this.  Return the copy.
-    const tvec3 operator-(const T& rhs) const
-    {
-        return tvec3(*this) -= rhs;
-    }
 
-    // Component-wise subtraction of another vector from this.
-    // Return a reference to this.
-    tvec3& operator-=(const tvec3& rhs)
-    {
-        x_ -= rhs.x_;
-        y_ -= rhs.y_;
-        z_ -= rhs.z_;
-        return *this;
-    }
+    template<scalar T_RHS = T, size_t N_RHS = N, enum align A_RHS = A>
+    tvec<T,N,A>& operator*=(const tvec<T_RHS,N_RHS,A_RHS>& rhs) { (*this) = (*this) * rhs; return (*this); }
+    template<scalar T_RHS = T, size_t N_RHS = N, enum align A_RHS = A>
+    tvec<T,N,A>& operator/=(const tvec<T_RHS,N_RHS,A_RHS>& rhs) { (*this) = (*this) / rhs; return (*this); }
+    template<scalar T_RHS = T, size_t N_RHS = N, enum align A_RHS = A>
+    tvec<T,N,A>& operator+=(const tvec<T_RHS,N_RHS,A_RHS>& rhs) { (*this) = (*this) + rhs; return (*this); }
+    template<scalar T_RHS = T, size_t N_RHS = N, enum align A_RHS = A>
+    tvec<T,N,A>& operator-=(const tvec<T_RHS,N_RHS,A_RHS>& rhs) { (*this) = (*this) - rhs; return (*this); }
 
-    // Component-wise subtraction of another vector from a copy of this.
-    // Return the copy.
-    const tvec3 operator-(const tvec3& rhs) const
-    {
-        return tvec3(*this) -= rhs;
-    }
+    template<scalar T_RHS = T>
+    tvec<T,N,A>& operator*=(const T_RHS& rhs) { (*this) *= tvec<T,N,A>((T)rhs); return *this; }
+    template<scalar T_RHS = T>
+    tvec<T,N,A>& operator/=(const T_RHS& rhs) { (*this) /= tvec<T,N,A>((T)rhs); return *this; }
+    template<scalar T_RHS = T>
+    tvec<T,N,A>& operator+=(const T_RHS& rhs) { (*this) += tvec<T,N,A>((T)rhs); return *this; }
+    template<scalar T_RHS = T>
+    tvec<T,N,A>& operator-=(const T_RHS& rhs) { (*this) -= tvec<T,N,A>((T)rhs); return *this; }
 
     // Compute the length of this and return it.
     float length() const
@@ -411,261 +176,32 @@ class tvec3
     void normalize()
     {
         float l = length();
-        x_ /= l;
-        y_ /= l;
-        z_ /= l;
+	if(l != 0 && l != 1)
+		(*this) /= l;
     }
 
     // Compute the dot product of two vectors.
-    static T dot(const tvec3& v1, const tvec3& v2)
+    template<scalar T_V1 = float, scalar T_V2 = float, size_t N_V1 = 3, size_t N_V2 = 3, enum align A_V1 = align::adaptive, enum align A_V2 = align::adaptive, scalar T_DST = decltype((T_V1)1*(T_V2)1 + (T_V1)1*(T_V2)1)>
+    static T_DST dot(const tvec<T_V1,N_V1,A_V1>& v1, const tvec<T_V2,N_V2,A_V2>& v2)
     {
-        return (v1.x_ * v2.x_) + (v1.y_ * v2.y_) + (v1.z_ * v2.z_); 
+    	tvec<T_DST, std::min(N_V1,N_V2)> dst = v1 * v2;
+	return std::accumulate(dst.begin(), dst.end(), (T_DST)0);
     }
 
     // Compute the cross product of two vectors.
-    static tvec3 cross(const tvec3& u, const tvec3& v)
+    template<scalar T_U = float, scalar T_V = float, size_t N_U = 3, size_t N_V = 3, enum align A_U = align::adaptive, enum align A_V = align::adaptive, scalar T_DST = decltype((T_U)1*(T_V)1)>
+    static tvec<T_DST,3,A> cross(const tvec<T_U,N_U,A_U>& u, const tvec<T_V,N_V,A_V>& v)
     {
-        return tvec3((u.y_ * v.z_) - (u.z_ * v.y_),
-                    (u.z_ * v.x_) - (u.x_ * v.z_),
-                    (u.x_ * v.y_) - (u.y_ * v.x_));
+        return (u * v.yzx() - u.yzx() * v).yzx();
     }
-
-private:
-    T x_;
-    T y_;
-    T z_;
 };
 
-// A template class for creating, managing and operating on a 4-element vector
-// of any type you like (intended for built-in types, but as long as it 
-// supports the basic arithmetic and assignment operators, any type should
-// work).
-template<typename T>
-class tvec4
-{
-public:
-    tvec4() :
-        x_(0),
-        y_(0),
-        z_(0),
-        w_(0) {}
-    tvec4(const T t) :
-        x_(t),
-        y_(t),
-        z_(t),
-        w_(t) {}
-    tvec4(const T x, const T y, const T z, const T w) :
-        x_(x),
-        y_(y),
-        z_(z),
-        w_(w) {}
-    tvec4(const tvec4& v) :
-        x_(v.x_),
-        y_(v.y_),
-        z_(v.z_),
-        w_(v.w_) {}
-    ~tvec4() {}
-
-    // Print the elements of the vector to standard out.
-    // Really only useful for debug and test.
-    void print() const
-    {
-        std::cout << "| " << x_ << " " << y_ << " " << z_ << " " << w_ << " |" << std::endl;
-    }
-
-    // Allow raw data access for API calls and the like.
-    // For example, it is valid to pass a tvec4<float> into a call to
-    // the OpenGL command "glUniform4fv()".
-    operator const T*() const { return &x_;}
-
-    // Get and set access members for the individual elements.
-    const T x() const { return x_; }
-    const T y() const { return y_; }
-    const T z() const { return z_; }
-    const T w() const { return w_; }
-
-    void x(const T& val) { x_ = val; }
-    void y(const T& val) { y_ = val; }
-    void z(const T& val) { z_ = val; }
-    void w(const T& val) { w_ = val; }
-
-    // A direct assignment of 'rhs' to this.  Return a reference to this.
-    tvec4& operator=(const tvec4& rhs)
-    {
-        if (this != &rhs)
-        {
-            x_ = rhs.x_;
-            y_ = rhs.y_;
-            z_ = rhs.z_;
-            w_ = rhs.w_;
-        }
-        return *this;
-    }
-
-    // Divide this by a scalar.  Return a reference to this.
-    tvec4& operator/=(const T& rhs)
-    {
-        x_ /= rhs;
-        y_ /= rhs;
-        z_ /= rhs;
-        w_ /= rhs;
-        return *this;
-    }
-
-    // Divide a copy of this by a scalar.  Return the copy.
-    const tvec4 operator/(const T& rhs) const
-    {
-        return tvec4(*this) /= rhs;
-    }
-
-    // Component-wise divide of this by another vector.
-    // Return a reference to this.
-    tvec4& operator/=(const tvec4& rhs)
-    {
-        x_ /= rhs.x_;
-        y_ /= rhs.y_;
-        z_ /= rhs.z_;
-        w_ /= rhs.w_;
-        return *this;
-    }
-
-    // Component-wise divide of a copy of this by another vector.
-    // Return the copy.
-    const tvec4 operator/(const tvec4& rhs) const
-    {
-        return tvec4(*this) /= rhs;
-    }
-
-    // Multiply this by a scalar.  Return a reference to this.
-    tvec4& operator*=(const T& rhs)
-    {
-        x_ *= rhs;
-        y_ *= rhs;
-        z_ *= rhs;
-        w_ *= rhs;
-        return *this;
-    }
-
-    // Multiply a copy of this by a scalar.  Return the copy.
-    const tvec4 operator*(const T& rhs) const
-    {
-        return tvec4(*this) *= rhs;
-    }
-
-    // Component-wise multiply of this by another vector.
-    // Return a reference to this.
-    tvec4& operator*=(const tvec4& rhs)
-    {
-        x_ *= rhs.x_;
-        y_ *= rhs.y_;
-        z_ *= rhs.z_;
-        w_ *= rhs.w_;
-        return *this;
-    }
-
-    // Component-wise multiply of a copy of this by another vector.
-    // Return the copy.
-    const tvec4 operator*(const tvec4& rhs) const
-    {
-        return tvec4(*this) *= rhs;
-    }
-
-    // Add a scalar to this.  Return a reference to this.
-    tvec4& operator+=(const T& rhs)
-    {
-        x_ += rhs;
-        y_ += rhs;
-        z_ += rhs;
-        w_ += rhs;
-        return *this;
-    }
-
-    // Add a scalar to a copy of this.  Return the copy.
-    const tvec4 operator+(const T& rhs) const
-    {
-        return tvec4(*this) += rhs;
-    }
-
-    // Component-wise addition of another vector to this.
-    // Return a reference to this.
-    tvec4& operator+=(const tvec4& rhs)
-    {
-        x_ += rhs.x_;
-        y_ += rhs.y_;
-        z_ += rhs.z_;
-        w_ += rhs.w_;
-        return *this;
-    }
-
-    // Component-wise addition of another vector to a copy of this.
-    // Return the copy.
-    const tvec4 operator+(const tvec4& rhs) const
-    {
-        return tvec4(*this) += rhs;
-    }
-
-    // Subtract a scalar from this.  Return a reference to this.
-    tvec4& operator-=(const T& rhs)
-    {
-        x_ -= rhs;
-        y_ -= rhs;
-        z_ -= rhs;
-        w_ -= rhs;
-        return *this;
-    }
-
-    // Subtract a scalar from a copy of this.  Return the copy.
-    const tvec4 operator-(const T& rhs) const
-    {
-        return tvec4(*this) -= rhs;
-    }
-
-    // Component-wise subtraction of another vector from this.
-    // Return a reference to this.
-    tvec4& operator-=(const tvec4& rhs)
-    {
-        x_ -= rhs.x_;
-        y_ -= rhs.y_;
-        z_ -= rhs.z_;
-        w_ -= rhs.w_;
-        return *this;
-    }
-
-    // Component-wise subtraction of another vector from a copy of this.
-    // Return the copy.
-    const tvec4 operator-(const tvec4& rhs) const
-    {
-        return tvec4(*this) -= rhs;
-    }
-
-    // Compute the length of this and return it.
-    float length() const
-    {
-        return sqrt(dot(*this, *this));
-    }
-
-    // Make this a unit vector.
-    void normalize()
-    {
-        float l = length();
-        x_ /= l;
-        y_ /= l;
-        z_ /= l;
-        w_ /= l;
-    }
-
-    // Compute the dot product of two vectors.
-    static T dot(const tvec4& v1, const tvec4& v2)
-    {
-        return (v1.x_ * v2.x_) + (v1.y_ * v2.y_) + (v1.z_ * v2.z_) + (v1.w_ * v2.w_); 
-    }
-
-private:
-    T x_;
-    T y_;
-    T z_;
-    T w_;
-};
+template<scalar T, enum align A = align::adaptive>
+using tvec2 = tvec<T, 2, A>;
+template<scalar T, enum align A = align::adaptive>
+using tvec3 = tvec<T, 3, A>;
+template<scalar T, enum align A = align::adaptive>
+using tvec4 = tvec<T, 4, A>;
 
 //
 // Convenience typedefs.  These are here to present a homogeneous view of these
@@ -695,20 +231,8 @@ typedef tvec4<bool> bvec4;
 
 // Global operators to allow for things like defining a new vector in terms of
 // a product of a scalar and a vector
-template<typename T>
-const LibMatrix::tvec2<T> operator*(const T t, const LibMatrix::tvec2<T>& v)
-{
-    return v * t;
-}
-
-template<typename T>
-const LibMatrix::tvec3<T> operator*(const T t, const LibMatrix::tvec3<T>& v)
-{
-    return v * t;
-}
-
-template<typename T>
-const LibMatrix::tvec4<T> operator*(const T t, const LibMatrix::tvec4<T>& v)
+template<scalar T, size_t N>
+const LibMatrix::tvec<T,N> operator*(const T t, const LibMatrix::tvec<T,N>& v)
 {
     return v * t;
 }