diff --git a/build/__pycache__/amsbuildlib4.cpython-312.pyc b/build/__pycache__/amsbuildlib4.cpython-312.pyc
new file mode 100644
index 0000000..f584104
Binary files /dev/null and b/build/__pycache__/amsbuildlib4.cpython-312.pyc differ
diff --git a/build_linux64/libamsmathutil25.linux64.a b/build_linux64/libamsmathutil25.linux64.a
index 30a862c..cd61525 100644
Binary files a/build_linux64/libamsmathutil25.linux64.a and b/build_linux64/libamsmathutil25.linux64.a differ
diff --git a/build_linux64/objstore/amscpptemplate25a_src2.o b/build_linux64/objstore/amscpptemplate25a_src2.o
index 96254c6..03a41ce 100644
Binary files a/build_linux64/objstore/amscpptemplate25a_src2.o and b/build_linux64/objstore/amscpptemplate25a_src2.o differ
diff --git a/build_linux64/objstore/amscpptemplate25a_template.o b/build_linux64/objstore/amscpptemplate25a_template.o
index 06d6ac8..1fc5d71 100644
Binary files a/build_linux64/objstore/amscpptemplate25a_template.o and b/build_linux64/objstore/amscpptemplate25a_template.o differ
diff --git a/build_linux64/objstore/amsmathtuil25_test1.o b/build_linux64/objstore/amsmathtuil25_test1.o
index 2fbeb00..2ab5e26 100644
Binary files a/build_linux64/objstore/amsmathtuil25_test1.o and b/build_linux64/objstore/amsmathtuil25_test1.o differ
diff --git a/build_linux64/objstore/amsmathutiil25_random.o b/build_linux64/objstore/amsmathutiil25_random.o
new file mode 100644
index 0000000..4a28951
Binary files /dev/null and b/build_linux64/objstore/amsmathutiil25_random.o differ
diff --git a/build_linux64/objstore/amsmathutil25_amsarray.o b/build_linux64/objstore/amsmathutil25_amsarray.o
index aefbb24..df8dc90 100644
Binary files a/build_linux64/objstore/amsmathutil25_amsarray.o and b/build_linux64/objstore/amsmathutil25_amsarray.o differ
diff --git a/build_linux64/objstore/amsmathutil25_amsarray_sort.o b/build_linux64/objstore/amsmathutil25_amsarray_sort.o
index 389530e..b876810 100644
Binary files a/build_linux64/objstore/amsmathutil25_amsarray_sort.o and b/build_linux64/objstore/amsmathutil25_amsarray_sort.o differ
diff --git a/build_linux64/objstore/amsmathutil25_mathfns1.o b/build_linux64/objstore/amsmathutil25_mathfns1.o
new file mode 100644
index 0000000..13bd9c9
Binary files /dev/null and b/build_linux64/objstore/amsmathutil25_mathfns1.o differ
diff --git a/build_linux64/objstore/amsmathutil25_vec2.o b/build_linux64/objstore/amsmathutil25_vec2.o
new file mode 100644
index 0000000..93f1c86
Binary files /dev/null and b/build_linux64/objstore/amsmathutil25_vec2.o differ
diff --git a/build_linux64/objstore/amsmathutil25_vec2f.o b/build_linux64/objstore/amsmathutil25_vec2f.o
new file mode 100644
index 0000000..57543d2
Binary files /dev/null and b/build_linux64/objstore/amsmathutil25_vec2f.o differ
diff --git a/build_linux64/objstore/amsmathutil25_vec3.o b/build_linux64/objstore/amsmathutil25_vec3.o
new file mode 100644
index 0000000..6cbc0a1
Binary files /dev/null and b/build_linux64/objstore/amsmathutil25_vec3.o differ
diff --git a/build_linux64/objstore/amsmathutil25_vec3f.o b/build_linux64/objstore/amsmathutil25_vec3f.o
new file mode 100644
index 0000000..3ca6459
Binary files /dev/null and b/build_linux64/objstore/amsmathutil25_vec3f.o differ
diff --git a/build_linux64/objstore/amsmathutil25_vec4.o b/build_linux64/objstore/amsmathutil25_vec4.o
new file mode 100644
index 0000000..65910fc
Binary files /dev/null and b/build_linux64/objstore/amsmathutil25_vec4.o differ
diff --git a/build_linux64/objstore/amsmathutil25_vec4f.o b/build_linux64/objstore/amsmathutil25_vec4f.o
new file mode 100644
index 0000000..6eb0c53
Binary files /dev/null and b/build_linux64/objstore/amsmathutil25_vec4f.o differ
diff --git a/build_linux64/objstore/amsmathutil2t_complex128.o b/build_linux64/objstore/amsmathutil2t_complex128.o
new file mode 100644
index 0000000..be57f9f
Binary files /dev/null and b/build_linux64/objstore/amsmathutil2t_complex128.o differ
diff --git a/build_linux64/objstore/amsmathutil2t_complex64.o b/build_linux64/objstore/amsmathutil2t_complex64.o
new file mode 100644
index 0000000..e2bbf48
Binary files /dev/null and b/build_linux64/objstore/amsmathutil2t_complex64.o differ
diff --git a/build_linux64/tests b/build_linux64/tests
index 1b1d8d6..aed5fd0 100644
Binary files a/build_linux64/tests and b/build_linux64/tests differ
diff --git a/include/amsmathutil25/amsmathutil25.hpp b/include/amsmathutil25/amsmathutil25.hpp
index c059954..493683d 100644
--- a/include/amsmathutil25/amsmathutil25.hpp
+++ b/include/amsmathutil25/amsmathutil25.hpp
@@ -9,6 +9,7 @@
 #include <vector>
 #include <thread>
 #include <functional>
+#include <mutex>
 
 namespace ams
 {
diff --git a/include/amsmathutil25/math/amsmathutil25_complex128.hpp b/include/amsmathutil25/math/amsmathutil25_complex128.hpp
index 54e9c78..f700ac8 100644
--- a/include/amsmathutil25/math/amsmathutil25_complex128.hpp
+++ b/include/amsmathutil25/math/amsmathutil25_complex128.hpp
@@ -6,6 +6,77 @@ namespace ams
 namespace cmp
 {
 
+    class complex
+    {
+        public:
+        double r,i;
+
+        complex();
+        complex(double _r, double _i);
+        complex(double _r);
+        complex(const complex &other);
+        complex& operator=(const complex& other);
+        complex& operator=(const double other);
+
+        friend complex operator-(const complex& z); //negation sign
+
+        complex operator+(const complex& z);
+        complex operator-(const complex& z);
+        complex operator*(const complex& z);
+        complex operator/(const complex& z);
+
+        complex operator+(const double & z);
+        complex operator-(const double & z);
+        complex operator*(const double & z);
+        complex operator/(const double & z);
+
+        double& operator[](int& ind);
+        const double& operator[](const int& ind) const;
+     
+        //comparison operators
+        bool operator==(const complex& z);
+        bool operator!=(const complex& z);
+        bool operator>(const complex& z);
+        bool operator<(const complex& z);
+        bool operator>=(const complex& z);
+        bool operator<=(const complex& z);
+        bool isnan();
+        bool isinf();
+        
+
+
+        bool isreal();
+        bool isimag();
+        bool iszero();
+        double arg();
+        double mag();
+        const complex conj();
+    };
+
+    double arg(ams::cmp::complex z);
+
+    double real(complex z);
+    double imag(complex z);
+    complex sin(complex z);
+    complex cos(complex z);
+    complex tan(complex z);
+    complex exp(complex z);
+    complex log(complex z);
+    double abs(complex z);
+    complex conj(complex z);
+
+    //need hyperbolic trig Functions
+    complex cosh(complex z);
+    complex sinh(complex z);
+    complex tanh(complex z);
+
+    complex pow(complex z1, complex z2);
+
+    //returns "complex sign" of complex number - 0, or a unit number with same argument
+    complex csgn(complex z);
+
+
+    typedef complex complex128;
     
 }; //end namespace cmp
 }; //end namespace ams
diff --git a/include/amsmathutil25/math/amsmathutil25_complex64.hpp b/include/amsmathutil25/math/amsmathutil25_complex64.hpp
index cd1f22a..40c275f 100644
--- a/include/amsmathutil25/math/amsmathutil25_complex64.hpp
+++ b/include/amsmathutil25/math/amsmathutil25_complex64.hpp
@@ -6,6 +6,74 @@ namespace ams
 namespace cmp
 {
 
+    class complex64
+    {
+        public:
+        float r,i;
+
+        complex64();
+        complex64(float _r, float _i);
+        complex64(float _r);
+        complex64(const complex64 &other);
+        complex64& operator=(const complex64& other);
+        complex64& operator=(const float other);
+
+        friend complex64 operator-(const complex64& z); //negation sign
+
+        complex64 operator+(const complex64& z);
+        complex64 operator-(const complex64& z);
+        complex64 operator*(const complex64& z);
+        complex64 operator/(const complex64& z);
+
+        complex64 operator+(const float & z);
+        complex64 operator-(const float & z);
+        complex64 operator*(const float & z);
+        complex64 operator/(const float & z);
+
+        float& operator[](int& ind);
+        const float& operator[](const int& ind) const;
+     
+        //comparison operators
+        bool operator==(const complex64& z);
+        bool operator!=(const complex64& z);
+        bool operator>(const complex64& z);
+        bool operator<(const complex64& z);
+        bool operator>=(const complex64& z);
+        bool operator<=(const complex64& z);
+        bool isnan();
+        bool isinf();
+        
+
+
+        bool isreal();
+        bool isimag();
+        bool iszero();
+        float arg();
+        float mag();
+        const complex64 conj();
+    };
+
+    float arg(ams::cmp::complex64 z);
+
+    float real(complex64 z);
+    float imag(complex64 z);
+    complex64 sin(complex64 z);
+    complex64 cos(complex64 z);
+    complex64 tan(complex64 z);
+    complex64 exp(complex64 z);
+    complex64 log(complex64 z);
+    float abs(complex64 z);
+    complex64 conj(complex64 z);
+
+    //need hyperbolic trig Functions
+    complex64 cosh(complex64 z);
+    complex64 sinh(complex64 z);
+    complex64 tanh(complex64 z);
+
+    complex64 pow(complex64 z1, complex64 z2);
+
+    //returns "complex64 sign" of complex64 number - 0, or a unit number with same argument
+    complex64 csgn(complex64 z);
     
 }; //end namespace cmp
 }; //end namespace ams
diff --git a/include/amsmathutil25/math/amsmathutil25_math.hpp b/include/amsmathutil25/math/amsmathutil25_math.hpp
index c64d798..bb20a19 100644
--- a/include/amsmathutil25/math/amsmathutil25_math.hpp
+++ b/include/amsmathutil25/math/amsmathutil25_math.hpp
@@ -4,6 +4,22 @@
 namespace ams
 {
 
+//OpenGL #defines a variable pi - the compiler doesn't like it if I #define a pi also...
+static const double pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164L;
+static const float pif = 3.1415926535897932384626433832795028841971693993751058209749445923078164;
+
+//this implementation doesn't work on all complilers?
+static const double nan = std::numeric_limits<double>::quiet_NaN();
+static const double inf = std::numeric_limits<double>::infinity();
+
+static const float fnan = std::numeric_limits<float>::quiet_NaN();
+static const float finf = std::numeric_limits<float>::infinity();
+
+static const double DLARGE = 1.0E300;
+static const double DSMALL = 1.0E-300;
+
+static const double vecmat_det_small = 1.0E-15;
+
 
 }; //end namespace ams
 
diff --git a/include/amsmathutil25/math/amsmathutil25_mathfns1.hpp b/include/amsmathutil25/math/amsmathutil25_mathfns1.hpp
index 27c6cb4..464f7b9 100644
--- a/include/amsmathutil25/math/amsmathutil25_mathfns1.hpp
+++ b/include/amsmathutil25/math/amsmathutil25_mathfns1.hpp
@@ -4,6 +4,71 @@
 namespace ams
 {
 
+bool isnan(double d);
+bool isinf(double d);
+bool isinfnan(double d);
+bool isnan(float f);
+bool isinf(float f);
+bool isinfnan(float f);
+
+double arg(double x, double y);
+float arg(float x, float y);
+
+double sgn(double d);
+float sgn(float f);
+
+
+//It has come to my attention that the behavior of Cs
+//fmod is odd: Negative numbers have negative moduli
+//I need to fix this with a proper modular function, whose values are in [0,N)
+double mod(double a, double b);
+float mod(float a, float b);
+int32_t mod(int32_t a, int32_t b);
+int64_t mod(int64_t a, int64_t b);
+
+int32_t truediv(int32_t x, int32_t y);
+int64_t truediv(int64_t x, int64_t y);
+
+double abs(double a);
+float abs(float a);
+int32_t abs(int32_t a);
+int64_t abs(int64_t a);
+
+template<typename T> T max(T a, T b);
+template<typename T> T min(T a, T b);
+template<typename T> T max(T a, T b, T c);
+template<typename T> T min(T a, T b, T c);
+
+//generic implementations
+template<typename T> T max(T a, T b)
+{
+    T ret = (a>b)? a : b;
+    return ret;
+}
+
+template<typename T> T min(T a, T b)
+{
+    T ret = (a<b)? a : b;
+    return ret;
+}
+
+template<typename T> T max(T a, T b, T c)
+{
+    T ret = a;
+    ret = (ret<b) ? b : ret;
+    ret = (ret<c) ? c : ret;
+    return ret;
+    
+}
+
+template<typename T> T min(T a, T b, T c)
+{
+    T ret = a;
+    ret = (ret>b) ? b : ret;
+    ret = (ret>c) ? c : ret;
+    return ret;
+}
+
 
 }; //end namespace ams
 
diff --git a/include/amsmathutil25/math/amsmathutil25_vec2.hpp b/include/amsmathutil25/math/amsmathutil25_vec2.hpp
index 189b42f..a6b9397 100644
--- a/include/amsmathutil25/math/amsmathutil25_vec2.hpp
+++ b/include/amsmathutil25/math/amsmathutil25_vec2.hpp
@@ -12,24 +12,90 @@ public:
 
     vec2();
     vec2(double _x, double _y);
-    vec2 operator=(vec2 rhs);
-
+    vec2(const vec2 &rhs);
+    vec2& operator=(const vec2 &rhs);
+    bool operator==(vec2 rhs);
+    
+    vec2 operator+(vec2 rhs) const;
+    vec2 operator-(vec2 rhs) const;
+    vec2 operator*(double rhs) const;
+    vec2 operator/(double rhs) const;
+    friend vec2 operator-(vec2 rhs);
 
+    double& operator[](int ind);
+    const double& operator[](int ind) const;
 };
 
+//vector operations
+double vec2_dot(vec2 v1, vec2 v2);
+double vec2_norm(vec2 v);
+vec2 vec2_normalize(vec2 v);
+double vec2_dist(vec2 v1, vec2 v2);
+double vec2_cross(vec2 v1, vec2 v2);
+vec2 vec2_proj(vec2 v1, vec2 v2);
+vec2 vec2_hodgedual(vec2 v); 
+double vec2_arg(vec2 v);
+
+//vector operations (old nomenclature)
+double vdot(vec2 v1, vec2 v2);
+double vnorm(vec2 v);
+vec2 vnormalize(vec2 v);
+double vdist(vec2 v1, vec2 v2);
+double vcross(vec2 v1, vec2 v2);
+vec2 vproj(vec2 v1, vec2 v2);
+
 class mat2
 {
 public:
     double data[4];
 
     mat2();
-    mat2(double _xx, double _xy, double _yx, double _yy);
+    mat2(double _xx, double _yx, double _xy, double _yy);
     mat2(const double *_data);
-    mat2(mat2 &rhs);
+    mat2(const mat2 &rhs);
+    mat2(const vec2 _col1, const vec2 _col2);
 
-    
+    mat2& operator=(const mat2 rhs);
+    bool operator==(const mat2 rhs);
+
+    mat2 operator+(mat2 rhs) const;
+    mat2 operator-(mat2 rhs) const;
+    mat2 operator*(mat2 rhs) const;
+    mat2 operator/(mat2 rhs) const;
+
+    mat2 operator*(double rhs) const;
+    mat2 operator/(double rhs) const;
+    vec2 operator*(vec2 rhs) const;
+
+    double& operator()(int I, int J);
+    double& operator[](int I);
+    const double& operator()(int I, int J) const;
+    const double& operator[](int I) const;
+    double& at(int I, int J);
+    const double& at(int I, int J) const;
+
+    double det();
+    mat2 inverse();
+    void invert();
+    mat2 transpose();
+    void dotranspose();
 };
 
+//matrix operations
+mat2 mat2_eye();
+mat2 mat2_zeros();
+mat2 mat2_ones();
+
+mat2 mat2_mult(mat2 a, mat2 b);
+vec2 mat2_mult(vec2 a, mat2 b);
+vec2 mat2_mult(mat2 a, vec2 b);
+
+double mat2_det(mat2 a);
+mat2 mat2_inverse(mat2 a);
+mat2 mat2_transpose(mat2 a);
+
+
+
 }; //end namespace ams
 
 #endif
diff --git a/include/amsmathutil25/math/amsmathutil25_vec2f.hpp b/include/amsmathutil25/math/amsmathutil25_vec2f.hpp
index 00ef248..646a058 100644
--- a/include/amsmathutil25/math/amsmathutil25_vec2f.hpp
+++ b/include/amsmathutil25/math/amsmathutil25_vec2f.hpp
@@ -4,6 +4,97 @@
 namespace ams
 {
 
+class vec2f
+{
+public:
+    float x;
+    float y;
+
+    vec2f();
+    vec2f(float _x, float _y);
+    vec2f(const vec2f &rhs);
+    vec2f& operator=(const vec2f &rhs);
+    bool operator==(vec2f rhs);
+    
+    vec2f operator+(vec2f rhs) const;
+    vec2f operator-(vec2f rhs) const;
+    vec2f operator*(float rhs) const;
+    vec2f operator/(float rhs) const;
+    friend vec2f operator-(vec2f rhs);
+
+    float& operator[](int ind);
+    const float& operator[](int ind) const;
+};
+
+//vector operations
+float vec2f_dot(vec2f v1, vec2f v2);
+float vec2f_norm(vec2f v);
+vec2f vec2f_normalize(vec2f v);
+float vec2f_dist(vec2f v1, vec2f v2);
+float vec2f_cross(vec2f v1, vec2f v2);
+vec2f vec2f_proj(vec2f v1, vec2f v2);
+vec2f vec2f_hodgedual(vec2f v); 
+float vec2f_arg(vec2f v);
+
+//vector operations (old nomenclature)
+float vdot(vec2f v1, vec2f v2);
+float vnorm(vec2f v);
+vec2f vnormalize(vec2f v);
+float vdist(vec2f v1, vec2f v2);
+float vcross(vec2f v1, vec2f v2);
+vec2f vproj(vec2f v1, vec2f v2);
+
+class mat2f
+{
+public:
+    float data[4];
+
+    mat2f();
+    mat2f(float _xx, float _yx, float _xy, float _yy);
+    mat2f(const float *_data);
+    mat2f(const mat2f &rhs);
+    mat2f(const vec2f _col1, const vec2f _col2);
+
+    mat2f& operator=(const mat2f rhs);
+    bool operator==(const mat2f rhs);
+
+    mat2f operator+(mat2f rhs) const;
+    mat2f operator-(mat2f rhs) const;
+    mat2f operator*(mat2f rhs) const;
+    mat2f operator/(mat2f rhs) const;
+
+    mat2f operator*(float rhs) const;
+    mat2f operator/(float rhs) const;
+    vec2f operator*(vec2f rhs) const;
+
+    float& operator()(int I, int J);
+    float& operator[](int I);
+    const float& operator()(int I, int J) const;
+    const float& operator[](int I) const;
+    float& at(int I, int J);
+    const float& at(int I, int J) const;
+
+    float det();
+    mat2f inverse();
+    void invert();
+    mat2f transpose();
+    void dotranspose();
+};
+
+//matrix operations
+mat2f mat2f_eye();
+mat2f mat2f_zeros();
+mat2f mat2f_ones();
+
+mat2f mat2f_mult(mat2f a, mat2f b);
+vec2f mat2f_mult(vec2f a, mat2f b);
+vec2f mat2f_mult(mat2f a, vec2f b);
+
+float mat2f_det(mat2f a);
+mat2f mat2f_inverse(mat2f a);
+mat2f mat2f_transpose(mat2f a);
+
+
 
 }; //end namespace ams
 
diff --git a/include/amsmathutil25/math/amsmathutil25_vec3.hpp b/include/amsmathutil25/math/amsmathutil25_vec3.hpp
index 3fca66e..7b06002 100644
--- a/include/amsmathutil25/math/amsmathutil25_vec3.hpp
+++ b/include/amsmathutil25/math/amsmathutil25_vec3.hpp
@@ -5,6 +5,99 @@ namespace ams
 {
 
 
+class vec3; 
+class mat3;
+
+class vec3
+{
+public:
+    double x;
+    double y;
+    double z;
+
+    vec3();
+    vec3(double _x, double _y, double _z);
+    vec3(const vec3 &rhs);
+    vec3& operator=(const vec3 &rhs);
+    bool operator==(vec3 rhs);
+    
+    vec3 operator+(vec3 rhs) const;
+    vec3 operator-(vec3 rhs) const;
+    vec3 operator*(double rhs) const;
+    vec3 operator/(double rhs) const;
+    friend vec3 operator-(vec3 rhs);
+
+    double& operator[](int ind);
+    const double& operator[](int ind) const;
+};
+
+//vector operations
+double vec3_dot(vec3 v1, vec3 v2);
+double vec3_norm(vec3 v);
+vec3 vec3_normalize(vec3 v);
+double vec3_dist(vec3 v1, vec3 v2);
+vec3 vec3_cross(vec3 v1, vec3 v2);
+vec3 vec3_proj(vec3 v1, vec3 v2);
+mat3 vec3_hodgedual(vec3 v); 
+
+//vector operations (old nomenclature)
+double vdot(vec3 v1, vec3 v2);
+double vnorm(vec3 v);
+vec3 vnormalize(vec3 v);
+double vdist(vec3 v1, vec3 v2);
+vec3 vcross(vec3 v1, vec3 v2);
+vec3 vproj(vec3 v1, vec3 v2);
+
+class mat3
+{
+public:
+    double data[9];
+
+    mat3();
+    mat3(double _xx, double _yx, double _zx, double _xy, double _yy, double _zy, double _xz, double _yz, double _zz);
+    mat3(const double *_data);
+    mat3(const mat3 &rhs);
+    mat3(const vec3 _col1, const vec3 _col2, const vec3 _col3);
+
+    mat3& operator=(const mat3 rhs);
+    bool operator==(const mat3 rhs);
+
+    mat3 operator+(mat3 rhs) const;
+    mat3 operator-(mat3 rhs) const;
+    mat3 operator*(mat3 rhs) const;
+    mat3 operator/(mat3 rhs) const;
+
+    mat3 operator*(double rhs) const;
+    mat3 operator/(double rhs) const;
+    vec3 operator*(vec3 rhs) const;
+
+    double& operator()(int I, int J);
+    double& operator[](int I);
+    const double& operator()(int I, int J) const;
+    const double& operator[](int I) const;
+    double& at(int I, int J);
+    const double& at(int I, int J) const;
+
+    double det();
+    mat3 inverse();
+    void invert();
+    mat3 transpose();
+    void dotranspose();
+};
+
+//matrix operations
+mat3 mat3_eye();
+mat3 mat3_zeros();
+mat3 mat3_ones();
+
+mat3 mat3_mult(mat3 a, mat3 b);
+vec3 mat3_mult(vec3 a, mat3 b);
+vec3 mat3_mult(mat3 a, vec3 b);
+
+double mat3_det(mat3 a);
+mat3 mat3_inverse(mat3 a);
+mat3 mat3_transpose(mat3 a);
+
 }; //end namespace ams
 
 #endif
diff --git a/include/amsmathutil25/math/amsmathutil25_vec3f.hpp b/include/amsmathutil25/math/amsmathutil25_vec3f.hpp
index 4d95651..16b6bfe 100644
--- a/include/amsmathutil25/math/amsmathutil25_vec3f.hpp
+++ b/include/amsmathutil25/math/amsmathutil25_vec3f.hpp
@@ -5,6 +5,99 @@ namespace ams
 {
 
 
+class vec3f; 
+class mat3f;
+
+class vec3f
+{
+public:
+    float x;
+    float y;
+    float z;
+
+    vec3f();
+    vec3f(float _x, float _y, float _z);
+    vec3f(const vec3f &rhs);
+    vec3f& operator=(const vec3f &rhs);
+    bool operator==(vec3f rhs);
+    
+    vec3f operator+(vec3f rhs) const;
+    vec3f operator-(vec3f rhs) const;
+    vec3f operator*(float rhs) const;
+    vec3f operator/(float rhs) const;
+    friend vec3f operator-(vec3f rhs);
+
+    float& operator[](int ind);
+    const float& operator[](int ind) const;
+};
+
+//vector operations
+float vec3f_dot(vec3f v1, vec3f v2);
+float vec3f_norm(vec3f v);
+vec3f vec3f_normalize(vec3f v);
+float vec3f_dist(vec3f v1, vec3f v2);
+vec3f vec3f_cross(vec3f v1, vec3f v2);
+vec3f vec3f_proj(vec3f v1, vec3f v2);
+mat3f vec3f_hodgedual(vec3f v); 
+
+//vector operations (old nomenclature)
+float vdot(vec3f v1, vec3f v2);
+float vnorm(vec3f v);
+vec3f vnormalize(vec3f v);
+float vdist(vec3f v1, vec3f v2);
+vec3f vcross(vec3f v1, vec3f v2);
+vec3f vproj(vec3f v1, vec3f v2);
+
+class mat3f
+{
+public:
+    float data[9];
+
+    mat3f();
+    mat3f(float _xx, float _yx, float _zx, float _xy, float _yy, float _zy, float _xz, float _yz, float _zz);
+    mat3f(const float *_data);
+    mat3f(const mat3f &rhs);
+    mat3f(const vec3f _col1, const vec3f _col2, const vec3f _col3);
+
+    mat3f& operator=(const mat3f rhs);
+    bool operator==(const mat3f rhs);
+
+    mat3f operator+(mat3f rhs) const;
+    mat3f operator-(mat3f rhs) const;
+    mat3f operator*(mat3f rhs) const;
+    mat3f operator/(mat3f rhs) const;
+
+    mat3f operator*(float rhs) const;
+    mat3f operator/(float rhs) const;
+    vec3f operator*(vec3f rhs) const;
+
+    float& operator()(int I, int J);
+    float& operator[](int I);
+    const float& operator()(int I, int J) const;
+    const float& operator[](int I) const;
+    float& at(int I, int J);
+    const float& at(int I, int J) const;
+
+    float det();
+    mat3f inverse();
+    void invert();
+    mat3f transpose();
+    void dotranspose();
+};
+
+//matrix operations
+mat3f mat3f_eye();
+mat3f mat3f_zeros();
+mat3f mat3f_ones();
+
+mat3f mat3f_mult(mat3f a, mat3f b);
+vec3f mat3f_mult(vec3f a, mat3f b);
+vec3f mat3f_mult(mat3f a, vec3f b);
+
+float mat3f_det(mat3f a);
+mat3f mat3f_inverse(mat3f a);
+mat3f mat3f_transpose(mat3f a);
+
 }; //end namespace ams
 
 #endif
diff --git a/include/amsmathutil25/math/amsmathutil25_vec4.hpp b/include/amsmathutil25/math/amsmathutil25_vec4.hpp
index 21264e1..daf90d3 100644
--- a/include/amsmathutil25/math/amsmathutil25_vec4.hpp
+++ b/include/amsmathutil25/math/amsmathutil25_vec4.hpp
@@ -4,6 +4,101 @@
 namespace ams
 {
 
+class vec4; 
+class mat4;
+
+class vec4
+{
+public:
+    double x;
+    double y;
+    double z;
+    double w;
+
+    vec4();
+    vec4(double _x, double _y, double _z, double _w);
+    vec4(const vec4 &rhs);
+    vec4& operator=(const vec4 &rhs);
+    bool operator==(vec4 rhs);
+    
+    vec4 operator+(vec4 rhs) const;
+    vec4 operator-(vec4 rhs) const;
+    vec4 operator*(double rhs) const;
+    vec4 operator/(double rhs) const;
+    friend vec4 operator-(vec4 rhs);
+
+    double& operator[](int ind);
+    const double& operator[](int ind) const;
+};
+
+//vector operations
+double vec4_dot(vec4 v1, vec4 v2);
+double vec4_norm(vec4 v);
+vec4 vec4_normalize(vec4 v);
+double vec4_dist(vec4 v1, vec4 v2);
+vec4 vec4_proj(vec4 v1, vec4 v2);
+
+//vector operations (old nomenclature)
+double vdot(vec4 v1, vec4 v2);
+double vnorm(vec4 v);
+vec4 vnormalize(vec4 v);
+double vdist(vec4 v1, vec4 v2);
+vec4 vproj(vec4 v1, vec4 v2);
+
+class mat4
+{
+public:
+    double data[16];
+
+    mat4();
+    mat4(
+        double _xx, double _yx, double _zx, double _wx,
+        double _xy, double _yy, double _zy, double _wy,
+        double _xz, double _yz, double _zz, double _wz,
+        double _xw, double _yw, double _zw, double _ww
+    );
+    mat4(const double *_data);
+    mat4(const mat4 &rhs);
+    mat4(const vec4 _col1, const vec4 _col2, const vec4 _col3, const vec4 _col4);
+
+    mat4& operator=(const mat4 rhs);
+    bool operator==(const mat4 rhs);
+
+    mat4 operator+(mat4 rhs) const;
+    mat4 operator-(mat4 rhs) const;
+    mat4 operator*(mat4 rhs) const;
+    mat4 operator/(mat4 rhs) const;
+
+    mat4 operator*(double rhs) const;
+    mat4 operator/(double rhs) const;
+    vec4 operator*(vec4 rhs) const;
+
+    double& operator()(int I, int J);
+    double& operator[](int I);
+    const double& operator()(int I, int J) const;
+    const double& operator[](int I) const;
+    double& at(int I, int J);
+    const double& at(int I, int J) const;
+
+    double det();
+    mat4 inverse();
+    void invert();
+    mat4 transpose();
+    void dotranspose();
+};
+
+//matrix operations
+mat4 mat4_eye();
+mat4 mat4_zeros();
+mat4 mat4_ones();
+
+mat4 mat4_mult(mat4 a, mat4 b);
+vec4 mat4_mult(vec4 a, mat4 b);
+vec4 mat4_mult(mat4 a, vec4 b);
+
+double mat4_det(mat4 a);
+mat4 mat4_inverse(mat4 a);
+mat4 mat4_transpose(mat4 a);
 
 }; //end namespace ams
 
diff --git a/include/amsmathutil25/math/amsmathutil25_vec4f.hpp b/include/amsmathutil25/math/amsmathutil25_vec4f.hpp
index 631840e..feda1c6 100644
--- a/include/amsmathutil25/math/amsmathutil25_vec4f.hpp
+++ b/include/amsmathutil25/math/amsmathutil25_vec4f.hpp
@@ -4,6 +4,101 @@
 namespace ams
 {
 
+class vec4f; 
+class mat4f;
+
+class vec4f
+{
+public:
+    float x;
+    float y;
+    float z;
+    float w;
+
+    vec4f();
+    vec4f(float _x, float _y, float _z, float _w);
+    vec4f(const vec4f &rhs);
+    vec4f& operator=(const vec4f &rhs);
+    bool operator==(vec4f rhs);
+    
+    vec4f operator+(vec4f rhs) const;
+    vec4f operator-(vec4f rhs) const;
+    vec4f operator*(float rhs) const;
+    vec4f operator/(float rhs) const;
+    friend vec4f operator-(vec4f rhs);
+
+    float& operator[](int ind);
+    const float& operator[](int ind) const;
+};
+
+//vector operations
+float vec4f_dot(vec4f v1, vec4f v2);
+float vec4f_norm(vec4f v);
+vec4f vec4f_normalize(vec4f v);
+float vec4f_dist(vec4f v1, vec4f v2);
+vec4f vec4f_proj(vec4f v1, vec4f v2);
+
+//vector operations (old nomenclature)
+float vdot(vec4f v1, vec4f v2);
+float vnorm(vec4f v);
+vec4f vnormalize(vec4f v);
+float vdist(vec4f v1, vec4f v2);
+vec4f vproj(vec4f v1, vec4f v2);
+
+class mat4f
+{
+public:
+    float data[16];
+
+    mat4f();
+    mat4f(
+        float _xx, float _yx, float _zx, float _wx,
+        float _xy, float _yy, float _zy, float _wy,
+        float _xz, float _yz, float _zz, float _wz,
+        float _xw, float _yw, float _zw, float _ww
+    );
+    mat4f(const float *_data);
+    mat4f(const mat4f &rhs);
+    mat4f(const vec4f _col1, const vec4f _col2, const vec4f _col3, const vec4f _col4);
+
+    mat4f& operator=(const mat4f rhs);
+    bool operator==(const mat4f rhs);
+
+    mat4f operator+(mat4f rhs) const;
+    mat4f operator-(mat4f rhs) const;
+    mat4f operator*(mat4f rhs) const;
+    mat4f operator/(mat4f rhs) const;
+
+    mat4f operator*(float rhs) const;
+    mat4f operator/(float rhs) const;
+    vec4f operator*(vec4f rhs) const;
+
+    float& operator()(int I, int J);
+    float& operator[](int I);
+    const float& operator()(int I, int J) const;
+    const float& operator[](int I) const;
+    float& at(int I, int J);
+    const float& at(int I, int J) const;
+
+    float det();
+    mat4f inverse();
+    void invert();
+    mat4f transpose();
+    void dotranspose();
+};
+
+//matrix operations
+mat4f mat4f_eye();
+mat4f mat4f_zeros();
+mat4f mat4f_ones();
+
+mat4f mat4f_mult(mat4f a, mat4f b);
+vec4f mat4f_mult(vec4f a, mat4f b);
+vec4f mat4f_mult(mat4f a, vec4f b);
+
+float mat4f_det(mat4f a);
+mat4f mat4f_inverse(mat4f a);
+mat4f mat4f_transpose(mat4f a);
 
 }; //end namespace ams
 
diff --git a/include/amsmathutil25/random/amsmathutil25_random.hpp b/include/amsmathutil25/random/amsmathutil25_random.hpp
index 9389861..24dddcb 100644
--- a/include/amsmathutil25/random/amsmathutil25_random.hpp
+++ b/include/amsmathutil25/random/amsmathutil25_random.hpp
@@ -6,6 +6,25 @@ namespace ams
 namespace rand
 {
 
+typedef ::int32_t amsmu_randt1;
+typedef ::int64_t amsmu_randt2;
+
+static const amsmu_randt1 dpr32_mod = ( ((amsmu_randt1)1) << ((amsmu_randt1)30) ) - (amsmu_randt1)1;
+static const amsmu_randt1 dpr32_mult1 = ( (amsmu_randt1) 1201633 );
+static const amsmu_randt1 dpr32_add1 = ( (amsmu_randt1) 293482 );
+
+static const amsmu_randt2 dpr64_mod = ( ((amsmu_randt2)1) << ((amsmu_randt2)62) ) - (amsmu_randt2)1;
+static const amsmu_randt2 dpr64_mult1 = ( (amsmu_randt2) 1201633L );
+static const amsmu_randt2 dpr64_add1 = ( (amsmu_randt2) 293482L );
+
+extern amsmu_randt1 dpr32_rseed;
+extern amsmu_randt2 dpr64_rseed;
+
+amsmu_randt1 dpr32_nextseed(amsmu_randt1 seed);
+amsmu_randt2 dpr64_nextseed(amsmu_randt2 seed);
+
+
+
 }; //end namespace rand
 }; //end namespace ams
 
diff --git a/include/amsmathutil25/testing/amsmathutil25_testing.hpp b/include/amsmathutil25/testing/amsmathutil25_testing.hpp
index 7ffaa92..af2311b 100644
--- a/include/amsmathutil25/testing/amsmathutil25_testing.hpp
+++ b/include/amsmathutil25/testing/amsmathutil25_testing.hpp
@@ -9,6 +9,9 @@ namespace amsmathutil25
     void test_amsarray1();
     void test_amsarray2();
 
+    void test_select();
+    void test_amsarray_sort();
+
 }; //end namespace amsmathutil25
 }; //end namespace ams
 
diff --git a/include/amsmathutil25/util/amsmathutil25_amsarray.hpp b/include/amsmathutil25/util/amsmathutil25_amsarray.hpp
index c3c89c3..2a5f11b 100644
--- a/include/amsmathutil25/util/amsmathutil25_amsarray.hpp
+++ b/include/amsmathutil25/util/amsmathutil25_amsarray.hpp
@@ -78,12 +78,15 @@ namespace ams
         amsarray<T> subarray(amsarray_size_t ind1, amsarray_size_t ind2) const;
 
         //returns an array that is {thisarray[inds[0]],thisarray[inds[1]],...}
-        amsarray<T> select(amsarray<amsarray_size_t> inds);
+        amsarray<T> select(const amsarray<amsarray_size_t> &inds);
         
         //returns an array of indices that is a permutation which will sort
         //this array in ascending order
         amsarray<amsarray_size_t> sort_permutation();
 
+        //sorts this array
+        int sort();
+
         //returns an array that is this array in reverse order
         amsarray<T> reverse();
 
diff --git a/include/amsmathutil25/util/amsmathutil25_amsarray_impl.hpp b/include/amsmathutil25/util/amsmathutil25_amsarray_impl.hpp
index 70fc738..8fd6ab2 100644
--- a/include/amsmathutil25/util/amsmathutil25_amsarray_impl.hpp
+++ b/include/amsmathutil25/util/amsmathutil25_amsarray_impl.hpp
@@ -122,21 +122,19 @@ namespace ams
         amsarray_size_t I;
         int res;
     
-        if(this!=&other)
+        res = this->resize(other.size());
+        if(res!=amsarray_success)
         {
-            res = this->resize(other.size());
-            if(res!=amsarray_success)
-            {
-                return *this;
-            }
-            else
+            return *this;
+        }
+        else
+        {
+            for(I=0;I<length && I<other.size();I++)
             {
-                for(I=0;I<length && I<other.size();I++)
-                {
-                    data[I] = other[I];
-                }
+                data[I] = other[I];
             }
         }
+
         return *this;
     }
     
@@ -565,14 +563,14 @@ template<typename T> bool amsarray<T>::operator!=(const amsarray<T>& other) cons
 template<typename T> int amsarray<T>::append(const T& val)
 {
     int ret = amsarray_success;
-    ret = insert(val,length);
+    ret = insert(length,val);
     return ret;
 }
 
 template<typename T> int amsarray<T>::prepend(const T& val)
 {
     int ret = amsarray_success;
-    ret = insert(val,0);
+    ret = insert(0,val);
     return ret;
 }
 
@@ -603,8 +601,8 @@ template<typename T> T amsarray<T>::pop_front()
 template<typename T> void amsarray_select_tf(
     int threadnum,
     int nthreads,
-    amsarray<T> *array,
-    amsarray<amsarray_size_t> *inds,
+    const amsarray<T> *array,
+    const amsarray<amsarray_size_t> *inds,
     amsarray<T> *ret
 )
 {
@@ -634,7 +632,7 @@ template<typename T> void amsarray_select_tf(
 }
 
 //returns an array that is {thisarray[inds[0]],thisarray[inds[1]],...}
-template<typename T> amsarray<T> amsarray<T>::select(amsarray<amsarray_size_t> inds)
+template<typename T> amsarray<T> amsarray<T>::select(const amsarray<amsarray_size_t> &inds)
 {
     int res;
     amsarray<T> ret;
@@ -667,7 +665,7 @@ template<typename T> amsarray<T> amsarray<T>::select(amsarray<amsarray_size_t> i
     else
     {
         threaded_execute(
-            &amsarray_select_tf,
+            &amsarray_select_tf<T>,
             inds.length,
             this,
             &inds,
@@ -681,7 +679,7 @@ template<typename T> amsarray<T> amsarray<T>::select(amsarray<amsarray_size_t> i
 template<typename T> void amsarray_reverse_tf(
     int threadnum,
     int nthreads,
-    amsarray<T> *array,
+    const amsarray<T> *array,
     amsarray<T> *ret
 )
 {
@@ -727,7 +725,7 @@ template<typename T> amsarray<T> amsarray<T>::reverse()
     else
     {
         threaded_execute(
-            &amsarray_reverse_tf,
+            &amsarray_reverse_tf<T>,
             length,
             this,
             &ret
@@ -737,7 +735,10 @@ template<typename T> amsarray<T> amsarray<T>::reverse()
     return ret;
 }
 
-
+template<typename T> void amsarray<T>::print(bool newline, int printstyle)
+{
+    //empty method - specialize for each type
+}
 
 
 
diff --git a/include/amsmathutil25/util/amsmathutil25_amsarray_sortimpl.hpp b/include/amsmathutil25/util/amsmathutil25_amsarray_sortimpl.hpp
index a5e68a9..dce9ec5 100644
--- a/include/amsmathutil25/util/amsmathutil25_amsarray_sortimpl.hpp
+++ b/include/amsmathutil25/util/amsmathutil25_amsarray_sortimpl.hpp
@@ -4,12 +4,13 @@
 namespace ams
 {
 
+void amsarray_permutation_swap(amsarray<amsarray_size_t> *permutation, amsarray_size_t I, amsarray_size_t J);
+
 
 template<typename T> int amsarray_quicksort_round(
-    amsarray<T> *array, 
-    amsarray<amsarray_size_t> *permarray, 
-    ams::pair<amsarray_size_t,amsarray_size_t> range,
-    amsarray_size_t *pivot_index,
+    amsarray<T> *array,                   //size N - array to sort
+    amsarray<amsarray_size_t> *permarray, //size N - permutation of sorting
+    ams::pair<amsarray_size_t,amsarray_size_t> range, //range over which to quicksort
     ams::pair<amsarray_size_t,amsarray_size_t> *leftrange,
     ams::pair<amsarray_size_t,amsarray_size_t> *rightrange
 )
@@ -19,13 +20,14 @@ template<typename T> int amsarray_quicksort_round(
     amsarray_size_t I,J,P;
     T v1,v2;
     amsarray_size_t tmp;
+    
+
 
     b1 = range.a < 0 || range.b < 0;
     b2 = (range.b - range.a) < 2;
-    if((b1 || b2)
+    if(b1 || b2)
     {
         //there is no more work to be done within this range
-        *pivot_index = -1;
         *leftrange = ams::pair<amsarray_size_t,amsarray_size_t>(-1,-1);
         *rightrange = ams::pair<amsarray_size_t,amsarray_size_t>(-1,-1);
         ret = -1;
@@ -35,16 +37,15 @@ template<typename T> int amsarray_quicksort_round(
     {
         //two element range - sort directly
         v1 = array->data[permarray->data[range.a]];
-        v2 = array->data[permarray->data[range.b]];
+        v2 = array->data[permarray->data[range.b-1]];
         if(v2<v1)
         {
             //swap permutation indices
             tmp = permarray->data[range.a];
-            permarray->data[range.a] = permarray->data[range.b];
-            permarray->data[range.b] = tmp;
+            permarray->data[range.a] = permarray->data[range.b-1];
+            permarray->data[range.b-1] = tmp;
         }
         //there is no more work to be done within this range
-        *pivot_index = -1;
         *leftrange = ams::pair<amsarray_size_t,amsarray_size_t>(-1,-1);
         *rightrange = ams::pair<amsarray_size_t,amsarray_size_t>(-1,-1);
         ret = -1;
@@ -52,19 +53,261 @@ template<typename T> int amsarray_quicksort_round(
     }
     else
     {
-        //perform quicksort iteration
+        //perform quicksort round
 
         //choose midpoint pivot
         P = (range.a + range.b)/2;
         P = (P<range.a) ? range.a : P;
         P = (P>=range.b) ? range.b-1 : P;
-
+        //swap pivot to end of range
+        amsarray_permutation_swap(permarray,P,range.b-1);
+        P = range.b-1;
         
+        J = range.a;
+        for(I=range.a;I<range.b-1;I++)
+        {
+            if(array->data[permarray->data[I]]<array->data[permarray->data[P]])
+            {
+                if(J!=I)
+                {
+                    amsarray_permutation_swap(permarray,I,J);
+                    J++;
+                }
+                else
+                {
+                    J++;
+                }
+            }
+        }
+
+        if(array->data[permarray->data[J]]<array->data[permarray->data[P]])
+        {
+            J++;
+            amsarray_permutation_swap(permarray,P,J);
+
+        }
+        else
+        {
+            amsarray_permutation_swap(permarray,P,J);
+        }
+        P = J;
+
+        if(P-range.a<=0)
+        {
+            leftrange->a = -1;
+            leftrange->b = -1;
+        }
+        else
+        {
+            leftrange->a = range.a;
+            leftrange->b = P;
+        }
+        if(range.b-(P+1)<=0)
+        {
+            rightrange->a = -1;
+            rightrange->b = -1;
+        }
+        else
+        {
+            rightrange->a = P+1;
+            rightrange->b = range.b;
+        }
+
+        ret = 1; //there is more work to do
+    }
+
+    return ret;
+}
+
+template<typename T> int amsarray_quicksort_unthreaded(
+    amsarray<T> *array,                   //size N - array to sort
+    amsarray<amsarray_size_t> *permarray //size N - permutation of sorting
+)
+{
+    int ret = amsarray_success;
+    int res;
+
+    amsarray<ams::pair<amsarray_size_t,amsarray_size_t>> ranges;
+    amsarray_size_t rangeptr = 0;
+    ams::pair<amsarray_size_t,amsarray_size_t> range,rangeleft,rangeright;
+
+    if(permarray->length!=array->length)
+    {
+        *permarray = permutation_identity(array->length);
+        if(permarray->length!=array->length)
+        {
+            ret = amsarray_failure;
+            return ret;
+        }
+    }
+
+    ranges.append(ams::pair<amsarray_size_t,amsarray_size_t>(0,array->length));
+    rangeptr = 0;
+    while(rangeptr<ranges.length)
+    {
+        range = ranges[rangeptr];
+        rangeptr++;
+        amsarray_quicksort_round(array,permarray,range,&rangeleft,&rangeright);
+        if(rangeleft.a>=0 && rangeleft.b>rangeleft.a)
+        {
+            ranges.append(rangeleft);
+        }
+        if(rangeright.a>=0 && rangeright.b>rangeright.a)
+        {
+            ranges.append(rangeright);
+        }
+    }
+
+    return ret;
+}
+
+template<typename T> void amsarray_quicksort_tf(
+    amsarray<T> *array,                   //size N - array to sort
+    amsarray<amsarray_size_t> *permarray, //size N - permutation of sorting
+    ams::pair<amsarray_size_t,amsarray_size_t> range,
+    amsarray<ams::pair<amsarray_size_t,amsarray_size_t>> *ranges,
+    std::mutex* threadlock
+)
+{
+    int res;
+
+    ams::pair<amsarray_size_t,amsarray_size_t> rangeleft,rangeright;
+
+    res = amsarray_quicksort_round(
+        array,
+        permarray,
+        range,
+        &rangeleft,
+        &rangeright
+    );
+
+    { //scope wrapper for std::lock_guard
+        std::lock_guard<std::mutex> lock(*threadlock);
+        //critical section
+        if(rangeleft.a>=0 && rangeleft.b>rangeleft.a)
+        {
+            ranges->append(rangeleft);
+        }
+        if(rangeright.a>=0 && rangeright.b>rangeright.a)
+        {
+            ranges->append(rangeright);
+        }
+    }
+    //end critical section (end of function)
+
+    return;
+}
+
+template<typename T> int amsarray_quicksort_threaded(
+    amsarray<T> *array,                   //size N - array to sort
+    amsarray<amsarray_size_t> *permarray //size N - permutation of sorting
+)
+{
+    int ret = amsarray_success;
+    int res;
+
+    amsarray<ams::pair<amsarray_size_t,amsarray_size_t>> ranges;
+    amsarray_size_t rangeptr = 0;
+    amsarray_size_t I = 0;
+    ams::pair<amsarray_size_t,amsarray_size_t> range;
+    amsarray<std::thread*> threads;
+    std::mutex threadlock;
+
+    int maxthreads = std::thread::hardware_concurrency();
+    maxthreads = (maxthreads < 1) ? 1 : maxthreads;
+    maxthreads = (maxthreads>amsmathutil25_maxthreads) ? amsmathutil25_maxthreads : maxthreads;
+
+    int nthreads;
+
+    if(permarray->length!=array->length)
+    {
+        *permarray = permutation_identity(array->length);
+        if(permarray->length!=array->length)
+        {
+            ret = amsarray_failure;
+            return ret;
+        }
+    }
+
+    threads.resize(maxthreads);
+    threads.setall(NULL);
+    rangeptr = 0;
+    ranges.append(ams::pair<amsarray_size_t,amsarray_size_t>(0,array->length));
 
+    while(rangeptr<ranges.length)
+    {
+        //spawn up to the maximum number of threads
+        nthreads = ranges.length-rangeptr;
+        nthreads = (nthreads>maxthreads) ? maxthreads : nthreads;
+
+        for(I=0;I<nthreads;I++)
+        {
+            range = ranges[rangeptr+I];
+            threads[I] = new(std::nothrow) std::thread(
+                amsarray_quicksort_tf<T>,
+                array,permarray,
+                range,
+                &ranges,
+                &threadlock
+            );
+        }
+        for(I=0;I<nthreads;I++)
+        {
+            if(threads[I]==NULL)
+            {
+                printf("amsarray_quicksort_threaded: error: thread %ld failed to spawn\n",I);
+                ret = amsarray_failure;
+            }
+        }
+        for(I=0;I<nthreads;I++)
+        {
+            if(threads[I]!=NULL)
+            {
+                threads[I]->join();
+                delete threads[I];
+                threads[I] = NULL;
+            }
+            rangeptr++;
+        }
     }
 
+    return ret;
+}
+
+template<typename T> int amsarray_quicksort(
+    amsarray<T> *array,                   //size N - array to sort
+    amsarray<amsarray_size_t> *permarray //size N - permutation of sorting
+)
+{
+    int ret = amsarray_success;
+    int res;
 
+    if(permarray->length!=array->length)
+    {
+        *permarray = permutation_identity(array->length);
+        if(permarray->length!=array->length)
+        {
+            ret = amsarray_failure;
+            return ret;
+        }
+    }
 
+    if(array->length<amsmathutil25_threadpsz)
+    {
+        //perform unthreaded quicksort
+        ret = amsarray_quicksort_unthreaded(
+            array,
+            permarray
+        );
+    }
+    else
+    {
+        //perform threaded quicksort
+        ret = amsarray_quicksort_threaded(
+            array,
+            permarray
+        );
+    }
 
     return ret;
 }
@@ -77,9 +320,37 @@ template<typename T> amsarray<amsarray_size_t> amsarray<T>::sort_permutation()
     int res;
     amsarray<amsarray_size_t> ret;
 
+    ret = permutation_identity(this->length);
+    if(ret.length!=this->length)
+    {
+        printf("sort_permutation: error - permutation array failed to allocate.\n");
+    }
+
+    res = amsarray_quicksort(
+        this,&ret
+    );
+
     return ret;
 }
 
+//returns an array of indices that is a permutation which will sort
+//this array in ascending order
+template<typename T> int amsarray<T>::sort()
+{
+    int ret = amsarray_success;
+    amsarray<amsarray_size_t> perm;
+    perm = sort_permutation();
+    if(perm.length==this->length)
+    {
+        *this = this->select(perm);
+    }
+    else
+    {
+        ret = amsarray_failure;
+    }
+
+    return ret;
+}
 
 
 
diff --git a/include/amsmathutil25/util/amsmathutil25_bufferops.hpp b/include/amsmathutil25/util/amsmathutil25_bufferops.hpp
index d272583..dab004c 100644
--- a/include/amsmathutil25/util/amsmathutil25_bufferops.hpp
+++ b/include/amsmathutil25/util/amsmathutil25_bufferops.hpp
@@ -11,6 +11,12 @@ namespace ams
     //threaded set - sets N elements of bufffer to val
     template<typename T> int buffer_set(T* buffer, buffer_size_t indstart, buffer_size_t indstop, const T val);
 
+    //threaded set - sets N elements of bufffer to val
+    template<typename T> int buffer_set_unthreaded(T* buffer, buffer_size_t N, const T val);
+
+    //threaded set - sets N elements of bufffer to val
+    template<typename T> int buffer_set_unthreaded(T* buffer, buffer_size_t indstart, buffer_size_t indstop, const T val);
+
     //threaded copy of bufferfrom[offsetfrom + I] to bufferto[offsetto+I] I = [0,N)
     template<typename T1, typename T2> int buffer_cast_copy(T1* bufferto, const T2* bufferfrom, 
         buffer_size_t offsetto, buffer_size_t offsetfrom, buffer_size_t N);
@@ -18,6 +24,13 @@ namespace ams
     //threaded copy of bufferfrom to bufferto
     template<typename T1, typename T2> int buffer_cast_copy(T1* bufferto, const T2* bufferfrom, buffer_size_t N);
 
+    //threaded copy of bufferfrom[offsetfrom + I] to bufferto[offsetto+I] I = [0,N)
+    template<typename T1, typename T2> int buffer_cast_copy_unthreaded(T1* bufferto, const T2* bufferfrom, 
+        buffer_size_t offsetto, buffer_size_t offsetfrom, buffer_size_t N);
+
+    //threaded copy of bufferfrom to bufferto
+    template<typename T1, typename T2> int buffer_cast_copy_unthreaded(T1* bufferto, const T2* bufferfrom, buffer_size_t N);
+
 }; //end namespace ams
 
 #include <amsmathutil25/util/amsmathutil25_bufferops_impl.hpp>
diff --git a/include/amsmathutil25/util/amsmathutil25_bufferops_impl.hpp b/include/amsmathutil25/util/amsmathutil25_bufferops_impl.hpp
index 4768370..b65348a 100644
--- a/include/amsmathutil25/util/amsmathutil25_bufferops_impl.hpp
+++ b/include/amsmathutil25/util/amsmathutil25_bufferops_impl.hpp
@@ -40,6 +40,7 @@ namespace ams
     template<typename T> int buffer_set(T* buffer, buffer_size_t indstart, buffer_size_t indstop, const T val)
     {
         int ret = amsmathutil25_success;
+        int res;
         buffer_size_t psize;
         buffer_size_t I;
 
@@ -59,15 +60,51 @@ namespace ams
         else
         {
 
-            ams::threaded_execute(
+            res = ams::threaded_execute(
                 &buffer_set_tf<T>,(int64_t) psize,
                 buffer,indstart,indstop,val
             );
+            ret = res;
 
         }
         return ret;
     }
 
+    //threaded set - sets N elements of bufffer to val
+    template<typename T> int buffer_set_unthreaded(T* buffer, buffer_size_t N, const T val)
+    {
+        int ret = amsmathutil25_success;
+        buffer_size_t I;
+
+        for(I=0;I<N;I++)
+        {
+            buffer[I] = val;
+        }
+
+        return ret;
+    }
+
+    //threaded set - sets N elements of bufffer to val
+    template<typename T> int buffer_set_unthreaded(T* buffer, buffer_size_t indstart, buffer_size_t indstop, const T val)
+    {
+        int ret = amsmathutil25_success;
+        buffer_size_t I;
+        buffer_size_t psize;
+
+        indstart = (indstart<0) ? 0 : indstart;
+        indstop = (indstop<0) ? 0 : indstop;
+        indstop = (indstop<indstart) ? indstart : indstop;
+        psize = indstop-indstart;
+
+        for(I=0;I<psize;I++)
+        {
+            buffer[I+indstart] = val;
+        }
+
+        return ret;
+    }
+
+
     //threaded set - sets N elements of bufffer to val
     template<typename T> int buffer_set(T* buffer, buffer_size_t N, const T val)
     {
@@ -128,11 +165,12 @@ namespace ams
         buffer_size_t offsetto, buffer_size_t offsetfrom, buffer_size_t N)
     {
         int ret = amsmathutil25_success;
-
-        ams::threaded_execute(
+        int res;
+        res = ams::threaded_execute(
             &buffer_cast_copy_tf<T1,T2>,(int64_t) N,
             bufferto,bufferfrom,offsetto,offsetfrom,N
         );
+        ret = res;
 
         return ret;
     }
@@ -155,6 +193,36 @@ namespace ams
         return ret;
     }
 
+    //threaded copy of bufferfrom[offsetfrom + I] to bufferto[offsetto+I] I = [0,N)
+    template<typename T1, typename T2> int buffer_cast_copy_unthreaded(T1* bufferto, const T2* bufferfrom, 
+        buffer_size_t offsetto, buffer_size_t offsetfrom, buffer_size_t N)
+    {
+        int ret = amsmathutil25_success;
+        buffer_size_t I;
+
+        for(I=0;I<N;I++)
+        {
+            bufferto[I+offsetto] = (T1) bufferfrom[I+offsetfrom];
+        }
+
+        return ret;
+    }
+
+    //threaded copy of bufferfrom to bufferto
+    template<typename T1, typename T2> int buffer_cast_copy_unthreaded(T1* bufferto, const T2* bufferfrom, buffer_size_t N)
+    {
+        int ret = amsmathutil25_success;
+        buffer_size_t I;
+
+        for(I=0;I<N;I++)
+        {
+            bufferto[I] = (T1) bufferfrom[I];
+        }
+
+        return ret;
+    }
+
+
 }; //end namespace ams
 
 #endif
diff --git a/src/amsmathutil25/math/amsmathutil25_mathfns1.cpp b/src/amsmathutil25/math/amsmathutil25_mathfns1.cpp
new file mode 100644
index 0000000..5930127
--- /dev/null
+++ b/src/amsmathutil25/math/amsmathutil25_mathfns1.cpp
@@ -0,0 +1,214 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+namespace amsmathutil25
+{
+
+bool isnan(double d)
+{
+    bool ret = 0;
+    if(d!=d)
+    {
+        ret = 1;
+    }
+    return ret;
+}
+
+bool isinf(double d)
+{
+    bool ret = 0;
+    ret = std::isinf(d);
+    return ret;
+}
+
+bool isinfnan(double d)
+{
+    bool ret = 0;
+    ret = ams::isinf(d) || ams::isnan(d);
+    return ret;
+}
+
+bool isnan(float f)
+{
+    bool ret = 0;
+    if(f!=f)
+    {
+        ret = 1;
+    }
+    return ret;
+}
+
+bool isinf(float f)
+{
+    bool ret = 0;
+    ret = std::isinf(f);
+    return ret;
+}
+
+bool isinfnan(float f)
+{
+    bool ret = 0;
+    ret = ams::isinf(f) || ams::isnan(f);
+    return ret;
+}
+
+double arg(double x, double y)
+{
+    double ret = 0.0;
+    double z = ::sqrt(x*x+y*y);
+
+    if(z>1.0E-10)
+    {
+        ret = ::atan2(y,x);
+        if(ret<0.0) ret = ret + 2.0*ams::pi;
+    }
+    return ret;
+}
+
+float arg(float x, float y)
+{
+    float ret = 0.0f;
+    float z = ::sqrtf(x*x+y*y);
+
+    if(z>1.0E-10)
+    {
+        ret = ::atan2f(y,x);
+        if(ret<0.0f) ret = ret + ams::pif;
+    }
+    return ret;
+}
+
+double sgn(double d)
+{
+    double ret = 0.0;
+    ret = (d>0.0) ? 1.0 : ret;
+    ret = (d<0.0) ? -1.0 : ret;
+
+    return ret;
+}
+
+float sgn(float f)
+{
+    float ret = 0.0f;
+    ret = (f>0.0f) ? 1.0f : ret;
+    ret = (f<0.0f) ? -1.0f : ret;
+
+    return ret;
+}
+
+//It has come to my attention that the behavior of Cs
+//fmod is odd: Negative numbers have negative moduli
+//I need to fix this with a proper modular function, whose values are in [0,N)
+double mod(double x, double n)
+{
+    x = ::fmod(x,n);
+  	if(x<0)
+  	{
+  		x = x + n;
+  	}
+  	return x;
+}
+float mod(float x, float n)
+{
+    x = ::fmodf(x,n);
+  	if(x<0)
+  	{
+  		x = x + n;
+  	}
+  	return x;
+}
+
+int32_t mod(int32_t x, int32_t n)
+{
+    x = x % n;
+    if(x<0)
+    {
+        x = x + n;
+    }
+    return x;
+}
+
+int64_t mod(int64_t x, int64_t n)
+{
+    x = x % n;
+    if(x<0)
+    {
+        x = x + n;
+    }
+    return x;
+}
+
+int32_t truediv(int32_t x, int32_t y)
+{
+    int32_t z = 0;
+    if(x>=0 && y>0)
+    {
+        z = x/y;
+    }
+    else if(x<0 && y>0)
+    {
+        z = -((-x)/y) - 1;
+    }
+    else if(x>=0 && y<0)
+    {
+        z = -(x/(-y)) - 1;
+    }
+    else if(x<0 && y<0)
+    {
+        z = ((-x)/(-y));
+    }
+
+    return z;
+}
+
+int64_t truediv(int64_t x, int64_t y)
+{
+    int64_t z = 0;
+    if(x>=0 && y>0)
+    {
+        z = x/y;
+    }
+    else if(x<0 && y>0)
+    {
+        z = -((-x)/y) - 1;
+    }
+    else if(x>=0 && y<0)
+    {
+        z = -(x/(-y)) - 1;
+    }
+    else if(x<0 && y<0)
+    {
+        z = ((-x)/(-y));
+    }
+
+    return z;
+}
+
+double abs(double a)
+{
+    double ret = ::fabs(a);
+    return ret;
+
+}
+
+float abs(float a)
+{
+    float ret = ::fabsf(a);
+    return ret;
+}
+
+int32_t abs(int32_t a)
+{
+    int32_t ret = (a<0) ? -a : a;
+    return ret;
+}
+
+int64_t abs(int64_t a)
+{
+    int64_t ret = (a<0) ? -a : a;
+    return ret;
+}
+
+};
+};
\ No newline at end of file
diff --git a/src/amsmathutil25/math/amsmathutil25_vec2.cpp b/src/amsmathutil25/math/amsmathutil25_vec2.cpp
new file mode 100644
index 0000000..dc29b80
--- /dev/null
+++ b/src/amsmathutil25/math/amsmathutil25_vec2.cpp
@@ -0,0 +1,466 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+
+
+vec2::vec2()
+{
+   x = 0.0;
+   y = 0.0; 
+}
+
+vec2::vec2(double _x, double _y)
+{
+    x = _x;
+    y = _y;
+}
+
+vec2& vec2::operator=(const vec2 &rhs)
+{
+    if(this!=&rhs)
+    {
+        this->x = rhs.x;
+        this->y = rhs.y;
+    }
+
+    return *this;
+}
+
+bool vec2::operator==(vec2 rhs)
+{
+    bool ret = 1;
+    ret = ret & (rhs.x == x);
+    ret = ret & (rhs.y == y);
+    return ret;
+}
+
+vec2::vec2(const vec2 &rhs)
+{
+    this->x = rhs.x;
+    this->y = rhs.y;
+    return;
+}
+
+vec2 vec2::operator+(vec2 rhs) const
+{
+    vec2 ret;
+    ret.x = this->x + rhs.x;
+    ret.y = this->y + rhs.y;
+    return ret;
+}
+
+vec2 vec2::operator-(vec2 rhs) const
+{
+    vec2 ret;
+    ret.x = this->x - rhs.x;
+    ret.y = this->y - rhs.y;
+    return ret;
+}
+
+vec2 vec2::operator*(double rhs) const
+{
+    vec2 ret;
+    ret.x = this->x*rhs;
+    ret.y = this->y*rhs;
+    return ret;
+}
+
+vec2 vec2::operator/(double rhs) const
+{
+    vec2 ret;
+    ret.x = this->x/rhs;
+    ret.y = this->y/rhs;
+    return ret;
+}
+
+vec2 operator-(vec2 rhs)
+{
+    vec2 ret;
+    ret.x = -rhs.x;
+    ret.y = -rhs.y;
+    return ret;
+}
+
+double& vec2::operator[](int ind)
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    return x;
+}
+
+const double& vec2::operator[](int ind) const
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    return x;
+}
+
+// vec2 vector operations
+double vec2_dot(vec2 v1, vec2 v2)
+{
+    double ret = v1.x*v2.x + v1.y*v2.y;
+    return ret;
+}
+
+double vec2_norm(vec2 v)
+{
+    double ret = ::sqrt(v.x*v.x+v.y*v.y);
+    return ret;
+}
+
+vec2 vec2_normalize(vec2 v)
+{
+    vec2 ret = vec2(0,0);
+    double nrm = vec2_norm(v);
+    if(nrm>0.0)
+    {
+        ret = v/nrm;
+    }
+    return ret;
+}
+
+double vec2_dist(vec2 v1, vec2 v2)
+{
+    return vec2_norm(v1-v2);
+}
+
+double vec2_cross(vec2 v1, vec2 v2)
+{
+    double ret = v1.x*v2.y - v1.y*v2.x;
+    return ret;
+}
+
+vec2 vec2_proj(vec2 v1, vec2 v2)
+{
+    vec2 ret;
+    vec2 v2n = vec2_normalize(v2);
+    ret = v2n*vec2_dot(v1,v2n);
+    return ret;
+}
+
+vec2 vec2_hodgedual(vec2 v)
+{
+    vec2 ret;
+    ret.x = -v.y;
+    ret.y = v.x;
+    return ret;
+}
+
+double vec2_arg(vec2 v)
+{
+    double ret;
+    ret = ::atan2(v.y,v.x);
+    if(ret<0.0)
+    {
+        ret = ret + 2.0*ams::pi;
+    }
+    return ret;
+}
+
+// vec2 vector operations (old nomenclature)
+double vdot(vec2 v1, vec2 v2)
+{
+    return vec2_dot(v1,v2);
+}
+
+double vnorm(vec2 v)
+{
+    return vec2_norm(v);
+}
+
+vec2 vnormalize(vec2 v)
+{
+    return vec2_normalize(v);
+}
+
+double vdist(vec2 v1, vec2 v2)
+{
+    return vec2_dist(v1,v2);
+}
+
+double vcross(vec2 v1, vec2 v2)
+{
+    return vec2_cross(v1,v2);
+}
+
+vec2 vproj(vec2 v1, vec2 v2)
+{
+    return vec2_proj(v1,v2);
+}
+
+// mat2 constructors and operators
+mat2::mat2()
+{
+    int I;
+    for(I=0;I<4;I++) data[I] = 0.0;
+    return;
+}
+
+mat2::mat2(double _xx, double _yx, double _xy, double _yy)
+{
+    data[0] = _xx;
+    data[1] = _yx;
+    data[2] = _xy;
+    data[3] = _yy;
+    return;
+}
+
+mat2::mat2(const double *_data)
+{
+    int I;
+    for(I=0;I<4;I++) data[I] = _data[I];
+}
+
+mat2::mat2(const mat2 &rhs)
+{
+    int I;
+    for(I=0;I<4;I++) data[I] = rhs.data[I];
+}
+
+mat2::mat2(const vec2 _col1, const vec2 _col2)
+{
+    data[0] = _col1[0];
+    data[1] = _col1[1];
+    data[2] = _col2[0];
+    data[3] = _col2[1];
+}
+
+mat2& mat2::operator=(const mat2 rhs)
+{
+    int I;
+    for(I=0;I<4;I++)
+    {
+        data[I] = rhs.data[I];
+    }
+    return *this;
+}
+
+bool mat2::operator==(const mat2 rhs)
+{
+    bool ret = 1;
+    int I;
+    for(I=0;I<4;I++)
+    {
+        if(data[I]!=rhs.data[I])
+        {
+            ret = 0;
+            break;
+        }
+    }
+
+    return ret;
+}
+
+mat2 mat2::operator+(mat2 rhs) const
+{
+    mat2 ret;
+    int I;
+    for(I=0;I<4;I++)
+    {
+        ret.data[I] = this->data[I] + rhs.data[I];
+    }
+    return ret;
+}
+
+mat2 mat2::operator-(mat2 rhs) const
+{
+    mat2 ret;
+    int I;
+    for(I=0;I<4;I++)
+    {
+        ret.data[I] = this->data[I] - rhs.data[I];
+    }
+    return ret;
+}
+
+mat2 mat2::operator*(mat2 rhs) const
+{
+    mat2 ret = mat2(0,0,0,0);
+    int I,J,K;
+    for(I=0;I<2;I++)
+    for(J=0;J<2;J++)
+    for(K=0;K<2;K++)
+        ret(I,K) += + (*this)(I,J)*rhs(J,K);
+    return ret;
+}
+
+mat2 mat2::operator/(mat2 rhs) const
+{
+    mat2 rhsi = rhs.inverse();
+    return (*this)*rhsi;
+}
+
+mat2 mat2::operator*(double rhs) const
+{
+    mat2 ret;
+    int I;
+    for(I=0;I<4;I++) ret.data[I] = this->data[I]*rhs;
+    return ret;
+}
+
+mat2 mat2::operator/(double rhs) const
+{
+    mat2 ret;
+    int I;
+    for(I=0;I<4;I++) ret.data[I] = this->data[I]/rhs;
+    return ret;
+}
+
+vec2 mat2::operator*(vec2 rhs) const
+{
+    vec2 ret = vec2(0,0);
+    int I,J;
+    for(I=0;I<2;I++)
+    for(J=0;J<2;J++)
+        ret[I] += (*this)(I,J)*rhs[J];
+    return ret;
+}
+
+double& mat2::operator()(int I, int J)
+{
+    if(I==0&&J==0) return data[0];
+    if(I==1&&J==0) return data[1];
+    if(I==0&&J==1) return data[2];
+    if(I==1&&J==1) return data[3];
+    return data[0];
+}
+
+double& mat2::operator[](int I)
+{
+    return data[I];
+}
+
+const double& mat2::operator()(int I, int J) const
+{
+    if(I==0&&J==0) return data[0];
+    if(I==1&&J==0) return data[1];
+    if(I==0&&J==1) return data[2];
+    if(I==1&&J==1) return data[3];
+    return data[0];
+}
+
+const double& mat2::operator[](int I) const
+{
+    return data[I];
+}
+
+double& mat2::at(int I, int J)
+{
+    return (*this)(I,J);
+}
+const double& mat2::at(int I, int J) const
+{
+    return (*this)(I,J);
+}
+
+double mat2::det()
+{
+    double ret = 0.0;
+    ret = ret + (*this)(0,0)*(*this)(1,1);
+    ret = ret - (*this)(1,0)*(*this)(0,1);
+    return ret;
+}
+
+mat2 mat2::inverse()
+{
+    mat2 ret = mat2_zeros();
+    double dt = this->det();
+
+    if(dt>=vecmat_det_small || dt<=-vecmat_det_small)
+    {
+        ret(0,0) = (*this)(1,1)/dt;
+        ret(0,1) = -(*this)(0,1)/dt;
+        ret(1,0) = -(*this)(1,0)/dt;
+        ret(1,1) = (*this)(0,0)/dt;
+    }
+    return ret;
+}
+
+void mat2::invert()
+{
+    mat2 q = this->inverse();
+    *this = q;
+    return;
+}
+
+mat2 mat2::transpose()
+{
+    mat2 ret;
+    int I,J;
+    for(I=0;I<2;I++)
+    for(J=0;J<2;J++)
+        ret(I,J) = this->at(J,I);
+    return ret;
+}
+
+void mat2::dotranspose()
+{
+    mat2 t = this->transpose();
+    (*this) = t;
+    return;
+}
+
+// mat2 matrix operations
+mat2 mat2_eye()
+{
+    mat2 ret;
+    ret(0,0) = 1.0;
+    ret(1,1) = 1.0;
+    ret(0,1) = 0.0;
+    ret(1,0) = 0.0;
+    return ret;
+}
+
+mat2 mat2_zeros()
+{
+    mat2 ret;
+    int I;
+    for(I=0;I<4;I++) ret.data[I] = 0.0;
+    return ret;
+}
+
+mat2 mat2_ones()
+{
+    mat2 ret;
+    int I;
+    for(I=0;I<4;I++) ret.data[I] = 1.0;
+    return ret;
+}
+
+mat2 mat2_mult(mat2 a, mat2 b)
+{
+    return a*b;
+}
+
+vec2 mat2_mult(vec2 a, mat2 b)
+{
+    vec2 ret = vec2(0,0);
+    int I,J;
+    for(I=0;I<2;I++)
+    for(J=0;J<2;J++)
+        ret[J] += a[I]*b(I,J);
+    return ret;
+}
+
+vec2 mat2_mult(mat2 a, vec2 b)
+{
+    return a*b;
+}
+
+double mat2_det(mat2 a)
+{
+    return a.det();
+}
+
+mat2 mat2_inverse(mat2 a)
+{
+    return a.inverse();
+}
+
+mat2 mat2_transpose(mat2 a)
+{
+    return a.transpose();
+}
+
+}; // end namespace ams
diff --git a/src/amsmathutil25/math/amsmathutil25_vec2f.cpp b/src/amsmathutil25/math/amsmathutil25_vec2f.cpp
new file mode 100644
index 0000000..2c782fe
--- /dev/null
+++ b/src/amsmathutil25/math/amsmathutil25_vec2f.cpp
@@ -0,0 +1,463 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+
+vec2f::vec2f()
+{
+   x = 0.0f;
+   y = 0.0f; 
+}
+
+vec2f::vec2f(float _x, float _y)
+{
+    x = _x;
+    y = _y;
+}
+
+vec2f& vec2f::operator=(const vec2f &rhs)
+{
+    if(this!=&rhs)
+    {
+        this->x = rhs.x;
+        this->y = rhs.y;
+    }
+
+    return *this;
+}
+
+bool vec2f::operator==(vec2f rhs)
+{
+    bool ret = 1;
+    ret = ret & (rhs.x == x);
+    ret = ret & (rhs.y == y);
+    return ret;
+}
+
+vec2f::vec2f(const vec2f &rhs)
+{
+    this->x = rhs.x;
+    this->y = rhs.y;
+    return;
+}
+
+vec2f vec2f::operator+(vec2f rhs) const
+{
+    vec2f ret;
+    ret.x = this->x + rhs.x;
+    ret.y = this->y + rhs.y;
+    return ret;
+}
+
+vec2f vec2f::operator-(vec2f rhs) const
+{
+    vec2f ret;
+    ret.x = this->x - rhs.x;
+    ret.y = this->y - rhs.y;
+    return ret;
+}
+
+vec2f vec2f::operator*(float rhs) const
+{
+    vec2f ret;
+    ret.x = this->x*rhs;
+    ret.y = this->y*rhs;
+    return ret;
+}
+
+vec2f vec2f::operator/(float rhs) const
+{
+    vec2f ret;
+    ret.x = this->x/rhs;
+    ret.y = this->y/rhs;
+    return ret;
+}
+
+vec2f operator-(vec2f rhs)
+{
+    vec2f ret;
+    ret.x = -rhs.x;
+    ret.y = -rhs.y;
+    return ret;
+}
+
+float& vec2f::operator[](int ind)
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    return x;
+}
+
+const float& vec2f::operator[](int ind) const
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    return x;
+}
+
+float vec2f_dot(vec2f v1, vec2f v2)
+{
+    float ret = v1.x*v2.x + v1.y*v2.y;
+    return ret;
+}
+
+float vec2f_norm(vec2f v)
+{
+    float ret = ::sqrtf(v.x*v.x+v.y*v.y);
+    return ret;
+}
+
+vec2f vec2f_normalize(vec2f v)
+{
+    vec2f ret = vec2f(0.0f,0.0f);
+    float nrm = vec2f_norm(v);
+    if(nrm>0.0f)
+    {
+        ret = v/nrm;
+    }
+    return ret;
+}
+
+float vec2f_dist(vec2f v1, vec2f v2)
+{
+    return vec2f_norm(v1-v2);
+}
+
+float vec2f_cross(vec2f v1, vec2f v2)
+{
+    float ret = v1.x*v2.y - v1.y*v2.x;
+    return ret;
+}
+
+vec2f vec2f_proj(vec2f v1, vec2f v2)
+{
+    vec2f ret;
+    vec2f v2n = vec2f_normalize(v2);
+    ret = v2n*vec2f_dot(v1,v2n);
+    return ret;
+}
+
+vec2f vec2f_hodgedual(vec2f v)
+{
+    vec2f ret;
+    ret.x = -v.y;
+    ret.y = v.x;
+    return ret;
+}
+
+float vec2f_arg(vec2f v)
+{
+    float ret;
+    ret = ::atan2f(v.y,v.x);
+    if(ret<0.0f)
+    {
+        ret = ret + 2.0f*ams::pif;
+    }
+    return ret;
+}
+
+// vec2f vector operations (old nomenclature)
+float vdot(vec2f v1, vec2f v2)
+{
+    return vec2f_dot(v1,v2);
+}
+
+float vnorm(vec2f v)
+{
+    return vec2f_norm(v);
+}
+
+vec2f vnormalize(vec2f v)
+{
+    return vec2f_normalize(v);
+}
+
+float vdist(vec2f v1, vec2f v2)
+{
+    return vec2f_dist(v1,v2);
+}
+
+float vcross(vec2f v1, vec2f v2)
+{
+    return vec2f_cross(v1,v2);
+}
+
+vec2f vproj(vec2f v1, vec2f v2)
+{
+    return vec2f_proj(v1,v2);
+}
+
+mat2f::mat2f()
+{
+    int I;
+    for(I=0;I<4;I++) data[I] = 0.0f;
+    return;
+}
+
+mat2f::mat2f(float _xx, float _yx, float _xy, float _yy)
+{
+    data[0] = _xx;
+    data[1] = _yx;
+    data[2] = _xy;
+    data[3] = _yy;
+    return;
+}
+
+mat2f::mat2f(const float *_data)
+{
+    int I;
+    for(I=0;I<4;I++) data[I] = _data[I];
+}
+
+mat2f::mat2f(const mat2f &rhs)
+{
+    int I;
+    for(I=0;I<4;I++) data[I] = rhs.data[I];
+}
+
+mat2f::mat2f(const vec2f _col1, const vec2f _col2)
+{
+    data[0] = _col1[0];
+    data[1] = _col1[1];
+    data[2] = _col2[0];
+    data[3] = _col2[1];
+}
+
+mat2f& mat2f::operator=(const mat2f rhs)
+{
+    int I;
+    for(I=0;I<4;I++)
+    {
+        data[I] = rhs.data[I];
+    }
+    return *this;
+}
+
+bool mat2f::operator==(const mat2f rhs)
+{
+    bool ret = 1;
+    int I;
+    for(I=0;I<4;I++)
+    {
+        if(data[I]!=rhs.data[I])
+        {
+            ret = 0;
+            break;
+        }
+    }
+
+    return ret;
+}
+
+mat2f mat2f::operator+(mat2f rhs) const
+{
+    mat2f ret;
+    int I;
+    for(I=0;I<4;I++)
+    {
+        ret.data[I] = this->data[I] + rhs.data[I];
+    }
+    return ret;
+}
+
+mat2f mat2f::operator-(mat2f rhs) const
+{
+    mat2f ret;
+    int I;
+    for(I=0;I<4;I++)
+    {
+        ret.data[I] = this->data[I] - rhs.data[I];
+    }
+    return ret;
+}
+
+mat2f mat2f::operator*(mat2f rhs) const
+{
+    mat2f ret = mat2f(0.0f,0.0f,0.0f,0.0f);
+    int I,J,K;
+    for(I=0;I<2;I++)
+    for(J=0;J<2;J++)
+    for(K=0;K<2;K++)
+        ret(I,K) += + (*this)(I,J)*rhs(J,K);
+    return ret;
+}
+
+mat2f mat2f::operator/(mat2f rhs) const
+{
+    mat2f rhsi = rhs.inverse();
+    return (*this)*rhsi;
+}
+
+mat2f mat2f::operator*(float rhs) const
+{
+    mat2f ret;
+    int I;
+    for(I=0;I<4;I++) ret.data[I] = this->data[I]*rhs;
+    return ret;
+}
+
+mat2f mat2f::operator/(float rhs) const
+{
+    mat2f ret;
+    int I;
+    for(I=0;I<4;I++) ret.data[I] = this->data[I]/rhs;
+    return ret;
+}
+
+vec2f mat2f::operator*(vec2f rhs) const
+{
+    vec2f ret = vec2f(0.0f,0.0f);
+    int I,J;
+    for(I=0;I<2;I++)
+    for(J=0;J<2;J++)
+        ret[I] += (*this)(I,J)*rhs[J];
+    return ret;
+}
+
+float& mat2f::operator()(int I, int J)
+{
+    if(I==0&&J==0) return data[0];
+    if(I==1&&J==0) return data[1];
+    if(I==0&&J==1) return data[2];
+    if(I==1&&J==1) return data[3];
+    return data[0];
+}
+
+float& mat2f::operator[](int I)
+{
+    return data[I];
+}
+
+const float& mat2f::operator()(int I, int J) const
+{
+    if(I==0&&J==0) return data[0];
+    if(I==1&&J==0) return data[1];
+    if(I==0&&J==1) return data[2];
+    if(I==1&&J==1) return data[3];
+    return data[0];
+}
+
+const float& mat2f::operator[](int I) const
+{
+    return data[I];
+}
+
+float& mat2f::at(int I, int J)
+{
+    return (*this)(I,J);
+}
+const float& mat2f::at(int I, int J) const
+{
+    return (*this)(I,J);
+}
+
+float mat2f::det()
+{
+    float ret = 0.0f;
+    ret = ret + (*this)(0,0)*(*this)(1,1);
+    ret = ret - (*this)(1,0)*(*this)(0,1);
+    return ret;
+}
+
+mat2f mat2f::inverse()
+{
+    mat2f ret = mat2f_zeros();
+    float dt = this->det();
+
+    if(dt>=vecmat_det_small || dt<=-vecmat_det_small)
+    {
+        ret(0,0) = (*this)(1,1)/dt;
+        ret(0,1) = -(*this)(0,1)/dt;
+        ret(1,0) = -(*this)(1,0)/dt;
+        ret(1,1) = (*this)(0,0)/dt;
+    }
+    return ret;
+}
+
+void mat2f::invert()
+{
+    mat2f q = this->inverse();
+    *this = q;
+    return;
+}
+
+mat2f mat2f::transpose()
+{
+    mat2f ret;
+    int I,J;
+    for(I=0;I<2;I++)
+    for(J=0;J<2;J++)
+        ret(I,J) = this->at(J,I);
+    return ret;
+}
+
+void mat2f::dotranspose()
+{
+    mat2f t = this->transpose();
+    (*this) = t;
+    return;
+}
+
+// mat2f matrix operations
+mat2f mat2f_eye()
+{
+    mat2f ret;
+    ret(0,0) = 1.0f;
+    ret(1,1) = 1.0f;
+    ret(0,1) = 0.0f;
+    ret(1,0) = 0.0f;
+    return ret;
+}
+
+mat2f mat2f_zeros()
+{
+    mat2f ret;
+    int I;
+    for(I=0;I<4;I++) ret.data[I] = 0.0f;
+    return ret;
+}
+
+mat2f mat2f_ones()
+{
+    mat2f ret;
+    int I;
+    for(I=0;I<4;I++) ret.data[I] = 1.0f;
+    return ret;
+}
+
+mat2f mat2f_mult(mat2f a, mat2f b)
+{
+    return a*b;
+}
+
+vec2f mat2f_mult(vec2f a, mat2f b)
+{
+    vec2f ret = vec2f(0.0f,0.0f);
+    int I,J;
+    for(I=0;I<2;I++)
+    for(J=0;J<2;J++)
+        ret[J] += a[I]*b(I,J);
+    return ret;
+}
+
+vec2f mat2f_mult(mat2f a, vec2f b)
+{
+    return a*b;
+}
+
+float mat2f_det(mat2f a)
+{
+    return a.det();
+}
+
+mat2f mat2f_inverse(mat2f a)
+{
+    return a.inverse();
+}
+
+mat2f mat2f_transpose(mat2f a)
+{
+    return a.transpose();
+}
+
+}; // end namespace ams
\ No newline at end of file
diff --git a/src/amsmathutil25/math/amsmathutil25_vec3.cpp b/src/amsmathutil25/math/amsmathutil25_vec3.cpp
new file mode 100644
index 0000000..213fe4e
--- /dev/null
+++ b/src/amsmathutil25/math/amsmathutil25_vec3.cpp
@@ -0,0 +1,506 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+
+vec3::vec3()
+{
+    x = 0.0;
+    y = 0.0;
+    z = 0.0;
+    return;
+}
+
+vec3::vec3(double _x, double _y, double _z)
+{
+    x = _x;
+    y = _y;
+    z = _z;
+    return;
+}
+
+vec3::vec3(const vec3 &rhs)
+{
+    x = rhs.x;
+    y = rhs.y;
+    z = rhs.z;
+    return;
+}
+
+vec3& vec3::operator=(const vec3 &rhs)
+{
+    this->x = rhs.x;
+    this->y = rhs.y;
+    this->z = rhs.z;
+    return *this;
+}
+
+bool vec3::operator==(vec3 rhs)
+{
+    return ((this->x == rhs.x) && (this->y == rhs.y) && (this->z == rhs.z));
+}
+
+// vec3 arithmetic operators
+vec3 vec3::operator+(vec3 rhs) const
+{
+    vec3 ret;
+    ret.x = this->x + rhs.x;
+    ret.y = this->y + rhs.y;
+    ret.z = this->z + rhs.z;
+    return ret;
+}
+
+vec3 vec3::operator-(vec3 rhs) const
+{
+    vec3 ret;
+    ret.x = this->x - rhs.x;
+    ret.y = this->y - rhs.y;
+    ret.z = this->z - rhs.z;
+    return ret;
+}
+
+vec3 vec3::operator*(double rhs) const
+{
+    vec3 ret;
+    ret.x = this->x*rhs;
+    ret.y = this->y*rhs;
+    ret.z = this->z*rhs;
+    return ret;
+}
+
+vec3 vec3::operator/(double rhs) const
+{
+    vec3 ret;
+    ret.x = this->x/rhs;
+    ret.y = this->y/rhs;
+    ret.z = this->z/rhs;
+    return ret;
+}
+
+vec3 operator-(vec3 rhs)
+{
+    vec3 ret;
+    ret.x = -rhs.x;
+    ret.y = -rhs.y;
+    ret.z = -rhs.z;
+    return ret;
+}
+
+// vec3 subscript operators
+double& vec3::operator[](int ind)
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    else if(ind==2) return z;
+    return x;
+}
+
+const double& vec3::operator[](int ind) const
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    else if(ind==2) return z;
+    return x;
+}
+
+// vector operations (new nomenclature)
+double vec3_dot(vec3 v1, vec3 v2)
+{
+    double ret = 0.0;
+    ret += v1.x*v2.x;
+    ret += v1.y*v2.y;
+    ret += v1.z*v2.z;
+    return ret;
+}
+
+double vec3_norm(vec3 v)
+{
+    double ret = ::sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
+    return ret;
+}
+
+vec3 vec3_normalize(vec3 v)
+{
+    vec3 ret = vec3(0,0,0);
+    double nrm = vec3_norm(v);
+    if(nrm>0.0)
+    {
+        ret = v/nrm;
+    }
+    return ret;
+}
+
+double vec3_dist(vec3 v1, vec3 v2)
+{
+    return vec3_norm(v1-v2);
+}
+
+vec3 vec3_cross(vec3 v1, vec3 v2)
+{
+    vec3 ret;
+    ret[0] = v1[1]*v2[2] - v1[2]*v2[1];
+    ret[1] = v1[2]*v2[0] - v1[0]*v2[2];
+    ret[2] = v1[0]*v2[1] - v1[1]*v2[0];
+    return ret;
+}
+
+vec3 vec3_proj(vec3 v1, vec3 v2)
+{
+    vec3 ret;
+    vec3 v2n = vec3_normalize(v2);
+    ret = v2n*vec3_dot(v1,v2n);
+    return ret;
+}
+
+mat3 vec3_hodgedual(vec3 v)
+{
+    mat3 ret;
+
+    ret(0,1) = v[2];
+    ret(1,2) = v[0];
+    ret(2,0) = v[1];
+
+    ret(1,0) = -v[2];
+    ret(2,1) = -v[0];
+    ret(0,2) = -v[1];
+
+    ret(0,0) = 0.0;
+    ret(1,1) = 0.0;
+    ret(2,2) = 0.0;
+
+    return ret;
+}
+
+// vector operations (old nomenclature)
+double vdot(vec3 v1, vec3 v2)
+{
+    return vec3_dot(v1,v2);
+}
+
+double vnorm(vec3 v)
+{
+    return vec3_norm(v);
+}
+
+vec3 vnormalize(vec3 v)
+{
+    return vec3_normalize(v);
+}
+
+double vdist(vec3 v1, vec3 v2)
+{
+    return vec3_dist(v1,v2);
+}
+
+vec3 vcross(vec3 v1, vec3 v2)
+{
+    return vec3_cross(v1,v2);
+}
+
+vec3 vproj(vec3 v1, vec3 v2)
+{
+    return vec3_proj(v1,v2);
+}
+
+// mat3 constructors
+mat3::mat3()
+{
+    int I;
+    for(I=0;I<9;I++) data[I] = 0.0;
+    return;
+}
+
+mat3::mat3(double _xx, double _yx, double _zx, double _xy, double _yy, double _zy, double _xz, double _yz, double _zz)
+{
+    data[0] = _xx;
+    data[1] = _yx;
+    data[2] = _zx;
+    
+    data[3] = _xy;
+    data[4] = _yy;
+    data[5] = _zy;
+    
+    data[6] = _xz;
+    data[7] = _yz;
+    data[8] = _zz;
+    
+    return;
+}
+
+mat3::mat3(const double *_data)
+{
+    int I;
+    for(I=0;I<9;I++) data[I] = _data[I];
+    return;
+}
+
+mat3::mat3(const mat3 &rhs)
+{
+    int I;
+    for(I=0;I<9;I++) data[I] = rhs.data[I];
+    return;
+}
+
+mat3::mat3(const vec3 _col1, const vec3 _col2, const vec3 _col3)
+{
+    data[0] = _col1[0];
+    data[1] = _col1[1];
+    data[2] = _col1[2];
+
+    data[3] = _col2[0];
+    data[4] = _col2[1];
+    data[5] = _col2[2];
+
+    data[6] = _col3[0];
+    data[7] = _col3[1];
+    data[8] = _col3[2];
+    
+    return;
+}
+
+// mat3 assignment and comparison operators
+mat3& mat3::operator=(const mat3 rhs)
+{
+    int I;
+    for(I=0;I<9;I++) data[I] = rhs.data[I];
+    return *this;
+}
+
+bool mat3::operator==(const mat3 rhs)
+{
+    bool ret = 1;
+    int I;
+    for(I=0;I<9;I++)
+    {
+        if(this->data[I]!=rhs.data[I])
+        {
+            ret = 0;
+            break;
+        }
+    }
+
+    return ret;
+}
+
+// mat3 arithmetic operators
+mat3 mat3::operator+(mat3 rhs) const
+{
+    mat3 ret;
+    int I;
+    for(I=0;I<9;I++)
+    {
+        ret.data[I] = this->data[I] + rhs.data[I];
+    }
+    return ret;
+}
+
+mat3 mat3::operator-(mat3 rhs) const
+{
+    mat3 ret;
+    int I;
+    for(I=0;I<9;I++)
+    {
+        ret.data[I] = this->data[I] - rhs.data[I];
+    }
+    return ret;
+}
+
+mat3 mat3::operator*(mat3 rhs) const
+{
+    mat3 ret = mat3_zeros();
+    int I,J,K;
+    for(I=0;I<3;I++)
+    for(J=0;J<3;J++)
+    for(K=0;K<3;K++)
+        ret(I,K) += + (*this)(I,J)*rhs(J,K);
+    return ret;
+}
+
+mat3 mat3::operator/(mat3 rhs) const
+{
+    mat3 rhsi = rhs.inverse();
+    return (*this)*rhsi;
+}
+
+mat3 mat3::operator*(double rhs) const
+{
+    mat3 ret;
+    int I;
+    for(I=0;I<9;I++) ret.data[I] = this->data[I]*rhs;
+    return ret;
+}
+
+mat3 mat3::operator/(double rhs) const
+{
+    mat3 ret;
+    int I;
+    for(I=0;I<9;I++) ret.data[I] = this->data[I]/rhs;
+    return ret;
+}
+
+vec3 mat3::operator*(vec3 rhs) const
+{
+    vec3 ret = vec3(0,0,0);
+    int I,J;
+    for(I=0;I<3;I++)
+    for(J=0;J<3;J++)
+        ret[I] += (*this)(I,J)*rhs[J];
+    return ret;
+}
+
+// mat3 access operators
+double& mat3::operator()(int I, int J)
+{
+    int ind = I + 3*J;
+    return data[ind];
+}
+
+double& mat3::operator[](int I)
+{
+    return data[I];
+}
+
+const double& mat3::operator()(int I, int J) const
+{
+    int ind = I + 3*J;
+    return data[ind];
+}
+
+const double& mat3::operator[](int I) const
+{
+    return data[I];
+}
+
+double& mat3::at(int I, int J)
+{
+    int ind = I + 3*J;
+    return data[ind];
+}
+
+const double& mat3::at(int I, int J) const
+{
+    int ind = I + 3*J;
+    return data[ind];
+}
+
+// mat3 member functions
+double mat3::det()
+{
+    double ret = 0.0;
+    ret += (*this)(0,0)*(*this)(1,1)*(*this)(2,2);
+    ret += (*this)(0,1)*(*this)(1,2)*(*this)(2,0);
+    ret += (*this)(0,2)*(*this)(1,0)*(*this)(2,1);
+
+    ret -= (*this)(2,0)*(*this)(1,1)*(*this)(0,2);
+    ret -= (*this)(2,1)*(*this)(1,2)*(*this)(0,0);
+    ret -= (*this)(2,2)*(*this)(1,0)*(*this)(0,1);
+    
+    return ret;
+}
+
+mat3 mat3::inverse()
+{
+    mat3 ret = mat3_zeros();
+	double dt = this->det();
+	if(dt>vecmat_det_small||dt<-vecmat_det_small)
+	{
+        ret(0,0) =  (this->at(1,1) * this->at(2,2) - this->at(1,2) * this->at(2,1))/dt;
+        ret(0,1) = -(this->at(1,0) * this->at(2,2) - this->at(1,2) * this->at(2,0))/dt;
+        ret(0,2) =  (this->at(1,0) * this->at(2,1) - this->at(1,1) * this->at(2,0))/dt;
+        ret(1,0) = -(this->at(0,1) * this->at(2,2) - this->at(0,2) * this->at(2,1))/dt;
+        ret(1,1) =  (this->at(0,0) * this->at(2,2) - this->at(0,2) * this->at(2,0))/dt;
+        ret(1,2) = -(this->at(0,0) * this->at(2,1) - this->at(0,1) * this->at(2,0))/dt;
+        ret(2,0) =  (this->at(0,1) * this->at(1,2) - this->at(0,2) * this->at(1,1))/dt;
+        ret(2,1) = -(this->at(0,0) * this->at(1,2) - this->at(0,2) * this->at(1,0))/dt;
+        ret(2,2) =  (this->at(0,0) * this->at(1,1) - this->at(0,1) * this->at(1,0))/dt;
+    }
+
+    return ret;
+}
+
+void mat3::invert()
+{
+    mat3 m = this->inverse();
+    *this = m;
+}
+
+mat3 mat3::transpose()
+{
+    mat3 ret;
+    int I,J;
+    for(I=0;I<3;I++)
+    for(J=0;J<3;J++)
+        ret(I,J) = this->at(J,I);
+    return ret;
+}
+
+void mat3::dotranspose()
+{
+    mat3 m = this->transpose();
+    *this = m;
+}
+
+// matrix operations (standalone functions)
+mat3 mat3_eye()
+{
+    mat3 ret = mat3_zeros();
+    ret(0,0) = 1.0;
+    ret(1,1) = 1.0;
+    ret(2,2) = 1.0;
+    return ret;
+}
+
+mat3 mat3_zeros()
+{
+    mat3 ret;
+    // for speed - the constructor already sets data[I] to 0.0
+    // int I;
+    // for(I=0;I<9;I++) ret.data[I] = 0.0;
+    return ret;
+}
+
+mat3 mat3_ones()
+{
+    mat3 ret;
+    int I;
+    for(I=0;I<9;I++) ret.data[I] = 1.0;
+    return ret;
+}
+
+mat3 mat3_mult(mat3 a, mat3 b)
+{
+    return a*b;
+}
+
+vec3 mat3_mult(vec3 a, mat3 b)
+{
+    vec3 ret = vec3(0,0,0);
+    int I,J;
+    for(I=0;I<3;I++)
+    for(J=0;J<3;J++)
+        ret[J] += a[I] * b(I,J);
+
+    return ret;
+}
+
+vec3 mat3_mult(mat3 a, vec3 b)
+{
+    return a*b;
+}
+
+double mat3_det(mat3 a)
+{
+    return a.det();
+}
+
+mat3 mat3_inverse(mat3 a)
+{
+    return a.inverse();
+}
+
+mat3 mat3_transpose(mat3 a)
+{
+    return a.transpose();
+}
+
+};
\ No newline at end of file
diff --git a/src/amsmathutil25/math/amsmathutil25_vec3f.cpp b/src/amsmathutil25/math/amsmathutil25_vec3f.cpp
new file mode 100644
index 0000000..75ab328
--- /dev/null
+++ b/src/amsmathutil25/math/amsmathutil25_vec3f.cpp
@@ -0,0 +1,506 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+
+vec3f::vec3f()
+{
+    x = 0.0;
+    y = 0.0;
+    z = 0.0;
+    return;
+}
+
+vec3f::vec3f(float _x, float _y, float _z)
+{
+    x = _x;
+    y = _y;
+    z = _z;
+    return;
+}
+
+vec3f::vec3f(const vec3f &rhs)
+{
+    x = rhs.x;
+    y = rhs.y;
+    z = rhs.z;
+    return;
+}
+
+vec3f& vec3f::operator=(const vec3f &rhs)
+{
+    this->x = rhs.x;
+    this->y = rhs.y;
+    this->z = rhs.z;
+    return *this;
+}
+
+bool vec3f::operator==(vec3f rhs)
+{
+    return ((this->x == rhs.x) && (this->y == rhs.y) && (this->z == rhs.z));
+}
+
+// vec3f arithmetic operators
+vec3f vec3f::operator+(vec3f rhs) const
+{
+    vec3f ret;
+    ret.x = this->x + rhs.x;
+    ret.y = this->y + rhs.y;
+    ret.z = this->z + rhs.z;
+    return ret;
+}
+
+vec3f vec3f::operator-(vec3f rhs) const
+{
+    vec3f ret;
+    ret.x = this->x - rhs.x;
+    ret.y = this->y - rhs.y;
+    ret.z = this->z - rhs.z;
+    return ret;
+}
+
+vec3f vec3f::operator*(float rhs) const
+{
+    vec3f ret;
+    ret.x = this->x*rhs;
+    ret.y = this->y*rhs;
+    ret.z = this->z*rhs;
+    return ret;
+}
+
+vec3f vec3f::operator/(float rhs) const
+{
+    vec3f ret;
+    ret.x = this->x/rhs;
+    ret.y = this->y/rhs;
+    ret.z = this->z/rhs;
+    return ret;
+}
+
+vec3f operator-(vec3f rhs)
+{
+    vec3f ret;
+    ret.x = -rhs.x;
+    ret.y = -rhs.y;
+    ret.z = -rhs.z;
+    return ret;
+}
+
+// vec3f subscript operators
+float& vec3f::operator[](int ind)
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    else if(ind==2) return z;
+    return x;
+}
+
+const float& vec3f::operator[](int ind) const
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    else if(ind==2) return z;
+    return x;
+}
+
+// vector operations (new nomenclature)
+float vec3f_dot(vec3f v1, vec3f v2)
+{
+    float ret = 0.0;
+    ret += v1.x*v2.x;
+    ret += v1.y*v2.y;
+    ret += v1.z*v2.z;
+    return ret;
+}
+
+float vec3f_norm(vec3f v)
+{
+    float ret = ::sqrtf(v.x*v.x+v.y*v.y+v.z*v.z);
+    return ret;
+}
+
+vec3f vec3f_normalize(vec3f v)
+{
+    vec3f ret = vec3f(0,0,0);
+    float nrm = vec3f_norm(v);
+    if(nrm>0.0)
+    {
+        ret = v/nrm;
+    }
+    return ret;
+}
+
+float vec3f_dist(vec3f v1, vec3f v2)
+{
+    return vec3f_norm(v1-v2);
+}
+
+vec3f vec3f_cross(vec3f v1, vec3f v2)
+{
+    vec3f ret;
+    ret[0] = v1[1]*v2[2] - v1[2]*v2[1];
+    ret[1] = v1[2]*v2[0] - v1[0]*v2[2];
+    ret[2] = v1[0]*v2[1] - v1[1]*v2[0];
+    return ret;
+}
+
+vec3f vec3f_proj(vec3f v1, vec3f v2)
+{
+    vec3f ret;
+    vec3f v2n = vec3f_normalize(v2);
+    ret = v2n*vec3f_dot(v1,v2n);
+    return ret;
+}
+
+mat3f vec3f_hodgedual(vec3f v)
+{
+    mat3f ret;
+
+    ret(0,1) = v[2];
+    ret(1,2) = v[0];
+    ret(2,0) = v[1];
+
+    ret(1,0) = -v[2];
+    ret(2,1) = -v[0];
+    ret(0,2) = -v[1];
+
+    ret(0,0) = 0.0;
+    ret(1,1) = 0.0;
+    ret(2,2) = 0.0;
+
+    return ret;
+}
+
+// vector operations (old nomenclature)
+float vdot(vec3f v1, vec3f v2)
+{
+    return vec3f_dot(v1,v2);
+}
+
+float vnorm(vec3f v)
+{
+    return vec3f_norm(v);
+}
+
+vec3f vnormalize(vec3f v)
+{
+    return vec3f_normalize(v);
+}
+
+float vdist(vec3f v1, vec3f v2)
+{
+    return vec3f_dist(v1,v2);
+}
+
+vec3f vcross(vec3f v1, vec3f v2)
+{
+    return vec3f_cross(v1,v2);
+}
+
+vec3f vproj(vec3f v1, vec3f v2)
+{
+    return vec3f_proj(v1,v2);
+}
+
+// mat3f constructors
+mat3f::mat3f()
+{
+    int I;
+    for(I=0;I<9;I++) data[I] = 0.0;
+    return;
+}
+
+mat3f::mat3f(float _xx, float _yx, float _zx, float _xy, float _yy, float _zy, float _xz, float _yz, float _zz)
+{
+    data[0] = _xx;
+    data[1] = _yx;
+    data[2] = _zx;
+    
+    data[3] = _xy;
+    data[4] = _yy;
+    data[5] = _zy;
+    
+    data[6] = _xz;
+    data[7] = _yz;
+    data[8] = _zz;
+    
+    return;
+}
+
+mat3f::mat3f(const float *_data)
+{
+    int I;
+    for(I=0;I<9;I++) data[I] = _data[I];
+    return;
+}
+
+mat3f::mat3f(const mat3f &rhs)
+{
+    int I;
+    for(I=0;I<9;I++) data[I] = rhs.data[I];
+    return;
+}
+
+mat3f::mat3f(const vec3f _col1, const vec3f _col2, const vec3f _col3)
+{
+    data[0] = _col1[0];
+    data[1] = _col1[1];
+    data[2] = _col1[2];
+
+    data[3] = _col2[0];
+    data[4] = _col2[1];
+    data[5] = _col2[2];
+
+    data[6] = _col3[0];
+    data[7] = _col3[1];
+    data[8] = _col3[2];
+    
+    return;
+}
+
+// mat3f assignment and comparison operators
+mat3f& mat3f::operator=(const mat3f rhs)
+{
+    int I;
+    for(I=0;I<9;I++) data[I] = rhs.data[I];
+    return *this;
+}
+
+bool mat3f::operator==(const mat3f rhs)
+{
+    bool ret = 1;
+    int I;
+    for(I=0;I<9;I++)
+    {
+        if(this->data[I]!=rhs.data[I])
+        {
+            ret = 0;
+            break;
+        }
+    }
+
+    return ret;
+}
+
+// mat3f arithmetic operators
+mat3f mat3f::operator+(mat3f rhs) const
+{
+    mat3f ret;
+    int I;
+    for(I=0;I<9;I++)
+    {
+        ret.data[I] = this->data[I] + rhs.data[I];
+    }
+    return ret;
+}
+
+mat3f mat3f::operator-(mat3f rhs) const
+{
+    mat3f ret;
+    int I;
+    for(I=0;I<9;I++)
+    {
+        ret.data[I] = this->data[I] - rhs.data[I];
+    }
+    return ret;
+}
+
+mat3f mat3f::operator*(mat3f rhs) const
+{
+    mat3f ret = mat3f_zeros();
+    int I,J,K;
+    for(I=0;I<3;I++)
+    for(J=0;J<3;J++)
+    for(K=0;K<3;K++)
+        ret(I,K) += + (*this)(I,J)*rhs(J,K);
+    return ret;
+}
+
+mat3f mat3f::operator/(mat3f rhs) const
+{
+    mat3f rhsi = rhs.inverse();
+    return (*this)*rhsi;
+}
+
+mat3f mat3f::operator*(float rhs) const
+{
+    mat3f ret;
+    int I;
+    for(I=0;I<9;I++) ret.data[I] = this->data[I]*rhs;
+    return ret;
+}
+
+mat3f mat3f::operator/(float rhs) const
+{
+    mat3f ret;
+    int I;
+    for(I=0;I<9;I++) ret.data[I] = this->data[I]/rhs;
+    return ret;
+}
+
+vec3f mat3f::operator*(vec3f rhs) const
+{
+    vec3f ret = vec3f(0,0,0);
+    int I,J;
+    for(I=0;I<3;I++)
+    for(J=0;J<3;J++)
+        ret[I] += (*this)(I,J)*rhs[J];
+    return ret;
+}
+
+// mat3f access operators
+float& mat3f::operator()(int I, int J)
+{
+    int ind = I + 3*J;
+    return data[ind];
+}
+
+float& mat3f::operator[](int I)
+{
+    return data[I];
+}
+
+const float& mat3f::operator()(int I, int J) const
+{
+    int ind = I + 3*J;
+    return data[ind];
+}
+
+const float& mat3f::operator[](int I) const
+{
+    return data[I];
+}
+
+float& mat3f::at(int I, int J)
+{
+    int ind = I + 3*J;
+    return data[ind];
+}
+
+const float& mat3f::at(int I, int J) const
+{
+    int ind = I + 3*J;
+    return data[ind];
+}
+
+// mat3f member functions
+float mat3f::det()
+{
+    float ret = 0.0;
+    ret += (*this)(0,0)*(*this)(1,1)*(*this)(2,2);
+    ret += (*this)(0,1)*(*this)(1,2)*(*this)(2,0);
+    ret += (*this)(0,2)*(*this)(1,0)*(*this)(2,1);
+
+    ret -= (*this)(2,0)*(*this)(1,1)*(*this)(0,2);
+    ret -= (*this)(2,1)*(*this)(1,2)*(*this)(0,0);
+    ret -= (*this)(2,2)*(*this)(1,0)*(*this)(0,1);
+    
+    return ret;
+}
+
+mat3f mat3f::inverse()
+{
+    mat3f ret = mat3f_zeros();
+	float dt = this->det();
+	if(dt>vecmat_det_small||dt<-vecmat_det_small)
+	{
+        ret(0,0) =  (this->at(1,1) * this->at(2,2) - this->at(1,2) * this->at(2,1))/dt;
+        ret(0,1) = -(this->at(1,0) * this->at(2,2) - this->at(1,2) * this->at(2,0))/dt;
+        ret(0,2) =  (this->at(1,0) * this->at(2,1) - this->at(1,1) * this->at(2,0))/dt;
+        ret(1,0) = -(this->at(0,1) * this->at(2,2) - this->at(0,2) * this->at(2,1))/dt;
+        ret(1,1) =  (this->at(0,0) * this->at(2,2) - this->at(0,2) * this->at(2,0))/dt;
+        ret(1,2) = -(this->at(0,0) * this->at(2,1) - this->at(0,1) * this->at(2,0))/dt;
+        ret(2,0) =  (this->at(0,1) * this->at(1,2) - this->at(0,2) * this->at(1,1))/dt;
+        ret(2,1) = -(this->at(0,0) * this->at(1,2) - this->at(0,2) * this->at(1,0))/dt;
+        ret(2,2) =  (this->at(0,0) * this->at(1,1) - this->at(0,1) * this->at(1,0))/dt;
+    }
+
+    return ret;
+}
+
+void mat3f::invert()
+{
+    mat3f m = this->inverse();
+    *this = m;
+}
+
+mat3f mat3f::transpose()
+{
+    mat3f ret;
+    int I,J;
+    for(I=0;I<3;I++)
+    for(J=0;J<3;J++)
+        ret(I,J) = this->at(J,I);
+    return ret;
+}
+
+void mat3f::dotranspose()
+{
+    mat3f m = this->transpose();
+    *this = m;
+}
+
+// matrix operations (standalone functions)
+mat3f mat3f_eye()
+{
+    mat3f ret = mat3f_zeros();
+    ret(0,0) = 1.0;
+    ret(1,1) = 1.0;
+    ret(2,2) = 1.0;
+    return ret;
+}
+
+mat3f mat3f_zeros()
+{
+    mat3f ret;
+    // for speed - the constructor already sets data[I] to 0.0
+    // int I;
+    // for(I=0;I<9;I++) ret.data[I] = 0.0;
+    return ret;
+}
+
+mat3f mat3f_ones()
+{
+    mat3f ret;
+    int I;
+    for(I=0;I<9;I++) ret.data[I] = 1.0;
+    return ret;
+}
+
+mat3f mat3f_mult(mat3f a, mat3f b)
+{
+    return a*b;
+}
+
+vec3f mat3f_mult(vec3f a, mat3f b)
+{
+    vec3f ret = vec3f(0,0,0);
+    int I,J;
+    for(I=0;I<3;I++)
+    for(J=0;J<3;J++)
+        ret[J] += a[I] * b(I,J);
+
+    return ret;
+}
+
+vec3f mat3f_mult(mat3f a, vec3f b)
+{
+    return a*b;
+}
+
+float mat3f_det(mat3f a)
+{
+    return a.det();
+}
+
+mat3f mat3f_inverse(mat3f a)
+{
+    return a.inverse();
+}
+
+mat3f mat3f_transpose(mat3f a)
+{
+    return a.transpose();
+}
+
+};
\ No newline at end of file
diff --git a/src/amsmathutil25/math/amsmathutil25_vec4.cpp b/src/amsmathutil25/math/amsmathutil25_vec4.cpp
new file mode 100644
index 0000000..9643d91
--- /dev/null
+++ b/src/amsmathutil25/math/amsmathutil25_vec4.cpp
@@ -0,0 +1,615 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+
+// vec4 constructors
+vec4::vec4()
+{
+    x = 0.0;
+    y = 0.0;
+    z = 0.0;
+    w = 0.0;
+    return;
+}
+
+vec4::vec4(double _x, double _y, double _z, double _w)
+{
+    x = _x;
+    y = _y;
+    z = _z;
+    w = _w;  
+    return;
+}
+
+vec4::vec4(const vec4 &rhs)
+{
+    x = rhs.x;
+    y = rhs.y;
+    z = rhs.z;
+    w = rhs.w;
+    return;
+}
+
+// vec4 assignment and comparison operators
+vec4& vec4::operator=(const vec4 &rhs)
+{
+    this->x = rhs.x;
+    this->y = rhs.y;
+    this->z = rhs.z;
+    this->w = rhs.w;
+    return *this;
+}
+
+bool vec4::operator==(vec4 rhs)
+{
+    return ((this->x == rhs.x) && (this->y == rhs.y) && (this->z == rhs.z) && (this->w == rhs.w));
+}
+
+// vec4 arithmetic operators
+vec4 vec4::operator+(vec4 rhs) const
+{
+    vec4 ret;
+    ret.x = this->x + rhs.x;
+    ret.y = this->y + rhs.y;
+    ret.z = this->z + rhs.z;
+    ret.w = this->w + rhs.w;
+    return ret;
+}
+
+vec4 vec4::operator-(vec4 rhs) const
+{
+    vec4 ret;
+    ret.x = this->x - rhs.x;
+    ret.y = this->y - rhs.y;
+    ret.z = this->z - rhs.z;
+    ret.w = this->w - rhs.w;
+    return ret;
+}
+
+vec4 vec4::operator*(double rhs) const
+{
+    vec4 ret;
+    ret.x = this->x*rhs;
+    ret.y = this->y*rhs;
+    ret.z = this->z*rhs;
+    ret.w = this->w*rhs;
+    return ret;
+}
+
+vec4 vec4::operator/(double rhs) const
+{
+    vec4 ret;
+    ret.x = this->x/rhs;
+    ret.y = this->y/rhs;
+    ret.z = this->z/rhs;
+    ret.w = this->w/rhs;
+    return ret;
+}
+
+vec4 operator-(vec4 rhs)
+{
+    vec4 ret;
+    ret.x = -rhs.x;
+    ret.y = -rhs.y;
+    ret.z = -rhs.z;
+    ret.w = -rhs.w;
+    return ret;
+}
+
+// vec4 subscript operators
+double& vec4::operator[](int ind)
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    else if(ind==2) return z;
+    else if(ind==3) return w;
+    return x;
+}
+
+const double& vec4::operator[](int ind) const
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    else if(ind==2) return z;
+    else if(ind==3) return w;
+    return x;
+}
+
+// vector operations (new nomenclature)
+double vec4_dot(vec4 v1, vec4 v2)
+{
+    double ret = 0.0;
+    ret += v1.x*v2.x;
+    ret += v1.y*v2.y;
+    ret += v1.z*v2.z;
+    ret += v1.w*v2.w;
+    return ret;
+}
+
+double vec4_norm(vec4 v)
+{
+    double ret = ::sqrt(v.x*v.x+v.y*v.y+v.z*v.z+v.w*v.w);
+    return ret;
+}
+
+vec4 vec4_normalize(vec4 v)
+{
+    vec4 ret = vec4(0,0,0,0);
+    double nrm = vec4_norm(v);
+    if(nrm>0.0)
+    {
+        ret = v/nrm;
+    }
+    return ret;
+}
+
+double vec4_dist(vec4 v1, vec4 v2)
+{
+    return vec4_norm(v1-v2);
+}
+
+vec4 vec4_proj(vec4 v1, vec4 v2)
+{
+    vec4 ret;
+    vec4 v2n = vec4_normalize(v2);
+    ret = v2n*vec4_dot(v1,v2n);
+    return ret;
+}
+
+// vector operations (old nomenclature)
+double vdot(vec4 v1, vec4 v2)
+{
+    return vec4_dot(v1,v2);
+}
+
+double vnorm(vec4 v)
+{
+    return vec4_norm(v);
+}
+
+vec4 vnormalize(vec4 v)
+{
+    return vec4_normalize(v);
+}
+
+double vdist(vec4 v1, vec4 v2)
+{
+    return vec4_dist(v1,v2);
+}
+
+vec4 vproj(vec4 v1, vec4 v2)
+{
+    return vec4_proj(v1,v2);
+}
+
+// mat4 constructors
+mat4::mat4()
+{
+    int I;
+    for(I=0;I<16;I++) data[I] = 0.0;
+    return;
+}
+
+mat4::mat4(
+    double _xx, double _yx, double _zx, double _wx,
+    double _xy, double _yy, double _zy, double _wy,
+    double _xz, double _yz, double _zz, double _wz,
+    double _xw, double _yw, double _zw, double _ww
+)
+{
+    data[0] = _xx;
+    data[1] = _yx;
+    data[2] = _zx;
+    data[3] = _wx;
+    
+    data[4] = _xy;
+    data[5] = _yy;
+    data[6] = _zy;
+    data[7] = _wy;
+    
+    data[8] = _xz;
+    data[9] = _yz;
+    data[10] = _zz;
+    data[11] = _wz;
+    
+    data[12] = _xw;
+    data[13] = _yw;
+    data[14] = _zw;
+    data[15] = _ww;
+    
+    return;
+}
+
+mat4::mat4(const double *_data)
+{
+    int I;
+    for(I=0;I<16;I++) data[I] = _data[I];
+    return;
+}
+
+mat4::mat4(const mat4 &rhs)
+{
+    int I;
+    for(I=0;I<16;I++) data[I] = rhs.data[I];
+    return;
+}
+
+mat4::mat4(const vec4 _col1, const vec4 _col2, const vec4 _col3, const vec4 _col4)
+{
+    data[0] = _col1[0];
+    data[1] = _col1[1];
+    data[2] = _col1[2];
+    data[3] = _col1[3];
+
+    data[4] = _col2[0];
+    data[5] = _col2[1];
+    data[6] = _col2[2];
+    data[7] = _col2[3];
+
+    data[8] = _col3[0];
+    data[9] = _col3[1];
+    data[10] = _col3[2];
+    data[11] = _col3[3];
+
+    data[12] = _col4[0];
+    data[13] = _col4[1];
+    data[14] = _col4[2];
+    data[15] = _col4[3];
+    
+    return;
+}
+
+// mat4 assignment and comparison operators
+mat4& mat4::operator=(const mat4 rhs)
+{
+    int I;
+    for(I=0;I<16;I++) data[I] = rhs.data[I];
+    return *this;
+}
+
+bool mat4::operator==(const mat4 rhs)
+{
+    bool ret = 1;
+    int I;
+    for(I=0;I<16;I++)
+    {
+        if(this->data[I]!=rhs.data[I])
+        {
+            ret = 0;
+            break;
+        }
+    }
+
+    return ret;
+}
+
+// mat4 arithmetic operators
+mat4 mat4::operator+(mat4 rhs) const
+{
+    mat4 ret;
+    int I;
+    for(I=0;I<16;I++)
+    {
+        ret.data[I] = this->data[I] + rhs.data[I];
+    }
+    return ret;
+}
+
+mat4 mat4::operator-(mat4 rhs) const
+{
+    mat4 ret;
+    int I;
+    for(I=0;I<16;I++)
+    {
+        ret.data[I] = this->data[I] - rhs.data[I];
+    }
+    return ret;
+}
+
+mat4 mat4::operator*(mat4 rhs) const
+{
+    mat4 ret = mat4_zeros();
+    int I,J,K;
+    for(I=0;I<4;I++)
+    for(J=0;J<4;J++)
+    for(K=0;K<4;K++)
+        ret(I,K) += (*this)(I,J)*rhs(J,K);
+    return ret;
+}
+
+mat4 mat4::operator/(mat4 rhs) const
+{
+    mat4 rhsi = rhs.inverse();
+    return (*this)*rhsi;
+}
+
+mat4 mat4::operator*(double rhs) const
+{
+    mat4 ret;
+    int I;
+    for(I=0;I<16;I++) ret.data[I] = this->data[I]*rhs;
+    return ret;
+}
+
+mat4 mat4::operator/(double rhs) const
+{
+    mat4 ret;
+    int I;
+    for(I=0;I<16;I++) ret.data[I] = this->data[I]/rhs;
+    return ret;
+}
+
+vec4 mat4::operator*(vec4 rhs) const
+{
+    vec4 ret = vec4(0,0,0,0);
+    int I,J;
+    for(I=0;I<4;I++)
+    for(J=0;J<4;J++)
+        ret[I] += (*this)(I,J)*rhs[J];
+    return ret;
+}
+
+// mat4 access operators
+double& mat4::operator()(int I, int J)
+{
+    int ind = I + 4*J;
+    return data[ind];
+}
+
+double& mat4::operator[](int I)
+{
+    return data[I];
+}
+
+const double& mat4::operator()(int I, int J) const
+{
+    int ind = I + 4*J;
+    return data[ind];
+}
+
+const double& mat4::operator[](int I) const
+{
+    return data[I];
+}
+
+double& mat4::at(int I, int J)
+{
+    int ind = I + 4*J;
+    return data[ind];
+}
+
+const double& mat4::at(int I, int J) const
+{
+    int ind = I + 4*J;
+    return data[ind];
+}
+
+double mat4::det()
+{
+    double ret = 0.0;
+    double a00,a01,a02,a03;
+    double a10,a11,a12,a13;
+    double a20,a21,a22,a23;
+    double a30,a31,a32,a33;
+  
+    a00 = this->at(0,0);
+    a01 = this->at(0,1);
+    a02 = this->at(0,2);
+    a03 = this->at(0,3);
+    a10 = this->at(1,0);
+    a11 = this->at(1,1);
+    a12 = this->at(1,2);
+    a13 = this->at(1,3);
+    a20 = this->at(2,0);
+    a21 = this->at(2,1);
+    a22 = this->at(2,2);
+    a23 = this->at(2,3);
+    a30 = this->at(3,0);
+    a31 = this->at(3,1);
+    a32 = this->at(3,2);
+    a33 = this->at(3,3);
+  
+    ret = a03*a12*a21*a30 -
+          a02*a13*a21*a30 -
+          a03*a11*a22*a30 +
+          a01*a13*a22*a30 +
+          a02*a11*a23*a30 -
+          a01*a12*a23*a30 -
+          a03*a12*a20*a31 +
+          a02*a13*a20*a31 +
+          a03*a10*a22*a31 -
+          a00*a13*a22*a31 -
+          a02*a10*a23*a31 +
+          a00*a12*a23*a31 +
+          a03*a11*a20*a32 -
+          a01*a13*a20*a32 -
+          a03*a10*a21*a32 +
+          a00*a13*a21*a32 +
+          a01*a10*a23*a32 -
+          a00*a11*a23*a32 -
+          a02*a11*a20*a33 +
+          a01*a12*a20*a33 +
+          a02*a10*a21*a33 -
+          a00*a12*a21*a33 -
+          a01*a10*a22*a33 +
+          a00*a11*a22*a33;
+
+    return ret;
+}
+
+mat4 mat4::inverse()
+{
+    mat4 ret = mat4_zeros();
+    double dt = this->det();
+    
+    double a00,a01,a02,a03;
+    double a10,a11,a12,a13;
+    double a20,a21,a22,a23;
+    double a30,a31,a32,a33;
+  
+    double b00,b01,b02,b03;
+    double b10,b11,b12,b13;
+    double b20,b21,b22,b23;
+    double b30,b31,b32,b33;
+  
+    a00 = this->at(0,0);
+    a01 = this->at(0,1);
+    a02 = this->at(0,2);
+    a03 = this->at(0,3);
+    a10 = this->at(1,0);
+    a11 = this->at(1,1);
+    a12 = this->at(1,2);
+    a13 = this->at(1,3);
+    a20 = this->at(2,0);
+    a21 = this->at(2,1);
+    a22 = this->at(2,2);
+    a23 = this->at(2,3);
+    a30 = this->at(3,0);
+    a31 = this->at(3,1);
+    a32 = this->at(3,2);
+    a33 = this->at(3,3);
+  
+    if(ams::abs(dt)>vecmat_det_small)
+    {
+      b00 = -a13*a22*a31 + a12*a23*a31 + a13*a21*a32 - a11*a23*a32 - a12*a21*a33 + a11*a22*a33;
+      b01 =  a03*a22*a31 - a02*a23*a31 - a03*a21*a32 + a01*a23*a32 + a02*a21*a33 - a01*a22*a33;
+      b02 = -a03*a12*a31 + a02*a13*a31 + a03*a11*a32 - a01*a13*a32 - a02*a11*a33 + a01*a12*a33;
+      b03 =  a03*a12*a21 - a02*a13*a21 - a03*a11*a22 + a01*a13*a22 + a02*a11*a23 - a01*a12*a23;
+      b10 =  a13*a22*a30 - a12*a23*a30 - a13*a20*a32 + a10*a23*a32 + a12*a20*a33 - a10*a22*a33;
+      b11 = -a03*a22*a30 + a02*a23*a30 + a03*a20*a32 - a00*a23*a32 - a02*a20*a33 + a00*a22*a33;
+      b12 =  a03*a12*a30 - a02*a13*a30 - a03*a10*a32 + a00*a13*a32 + a02*a10*a33 - a00*a12*a33;
+      b13 = -a03*a12*a20 + a02*a13*a20 + a03*a10*a22 - a00*a13*a22 - a02*a10*a23 + a00*a12*a23;
+      b20 = -a13*a21*a30 + a11*a23*a30 + a13*a20*a31 - a10*a23*a31 - a11*a20*a33 + a10*a21*a33;
+      b21 =  a03*a21*a30 - a01*a23*a30 - a03*a20*a31 + a00*a23*a31 + a01*a20*a33 - a00*a21*a33;
+      b22 = -a03*a11*a30 + a01*a13*a30 + a03*a10*a31 - a00*a13*a31 - a01*a10*a33 + a00*a11*a33;
+      b23 =  a03*a11*a20 - a01*a13*a20 - a03*a10*a21 + a00*a13*a21 + a01*a10*a23 - a00*a11*a23;
+      b30 =  a12*a21*a30 - a11*a22*a30 - a12*a20*a31 + a10*a22*a31 + a11*a20*a32 - a10*a21*a32;
+      b31 = -a02*a21*a30 + a01*a22*a30 + a02*a20*a31 - a00*a22*a31 - a01*a20*a32 + a00*a21*a32;
+      b32 =  a02*a11*a30 - a01*a12*a30 - a02*a10*a31 + a00*a12*a31 + a01*a10*a32 - a00*a11*a32;
+      b33 = -a02*a11*a20 + a01*a12*a20 + a02*a10*a21 - a00*a12*a21 - a01*a10*a22 + a00*a11*a22;
+      b00 = b00/dt;
+      b01 = b01/dt;
+      b02 = b02/dt;
+      b03 = b03/dt;
+      b10 = b10/dt;
+      b11 = b11/dt;
+      b12 = b12/dt;
+      b13 = b13/dt;
+      b20 = b20/dt;
+      b21 = b21/dt;
+      b22 = b22/dt;
+      b23 = b23/dt;
+      b30 = b30/dt;
+      b31 = b31/dt;
+      b32 = b32/dt;
+      b33 = b33/dt;
+      ret(0,0) = b00;
+      ret(0,1) = b01;
+      ret(0,2) = b02;
+      ret(0,3) = b03;
+      ret(1,0) = b10;
+      ret(1,1) = b11;
+      ret(1,2) = b12;
+      ret(1,3) = b13;
+      ret(2,0) = b20;
+      ret(2,1) = b21;
+      ret(2,2) = b22;
+      ret(2,3) = b23;
+      ret(3,0) = b30;
+      ret(3,1) = b31;
+      ret(3,2) = b32;
+      ret(3,3) = b33;
+    }
+    //this was STUPID. Gaah. Computer algebra system saves the day? I'd be surprised if this didn't end up *slower* than gaussian elimination. Don't do this again!
+
+    return ret;
+}
+
+void mat4::invert()
+{
+    mat4 m = this->inverse();
+    *this = m;
+}
+
+mat4 mat4::transpose()
+{
+    mat4 ret;
+    int I,J;
+    for(I=0;I<4;I++)
+    for(J=0;J<4;J++)
+        ret(I,J) = this->at(J,I);
+    return ret;
+}
+
+void mat4::dotranspose()
+{
+    mat4 m = this->transpose();
+    *this = m;
+}
+
+// matrix operations (standalone functions)
+mat4 mat4_eye()
+{
+    mat4 ret = mat4_zeros();
+    ret(0,0) = 1.0;
+    ret(1,1) = 1.0;
+    ret(2,2) = 1.0;
+    ret(3,3) = 1.0;
+    return ret;
+}
+
+mat4 mat4_zeros()
+{
+    mat4 ret;
+    // for speed - the constructor already sets data[I] to 0.0
+    // int I;
+    // for(I=0;I<16;I++) ret.data[I] = 0.0;
+    return ret;
+}
+
+mat4 mat4_ones()
+{
+    mat4 ret;
+    int I;
+    for(I=0;I<16;I++) ret.data[I] = 1.0;
+    return ret;
+}
+
+mat4 mat4_mult(mat4 a, mat4 b)
+{
+    return a*b;
+}
+
+vec4 mat4_mult(vec4 a, mat4 b)
+{
+    vec4 ret = vec4(0,0,0,0);
+    int I,J;
+    for(I=0;I<4;I++)
+    for(J=0;J<4;J++)
+        ret[J] += a[I] * b(I,J);
+
+    return ret;
+}
+
+vec4 mat4_mult(mat4 a, vec4 b)
+{
+    return a*b;
+}
+
+double mat4_det(mat4 a)
+{
+    return a.det();
+}
+
+mat4 mat4_inverse(mat4 a)
+{
+    return a.inverse();
+}
+
+mat4 mat4_transpose(mat4 a)
+{
+    return a.transpose();
+}
+
+};
\ No newline at end of file
diff --git a/src/amsmathutil25/math/amsmathutil25_vec4f.cpp b/src/amsmathutil25/math/amsmathutil25_vec4f.cpp
new file mode 100644
index 0000000..2b3ceb1
--- /dev/null
+++ b/src/amsmathutil25/math/amsmathutil25_vec4f.cpp
@@ -0,0 +1,615 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+
+// vec4f constructors
+vec4f::vec4f()
+{
+    x = 0.0f;
+    y = 0.0f;
+    z = 0.0f;
+    w = 0.0f;
+    return;
+}
+
+vec4f::vec4f(float _x, float _y, float _z, float _w)
+{
+    x = _x;
+    y = _y;
+    z = _z;
+    w = _w;  
+    return;
+}
+
+vec4f::vec4f(const vec4f &rhs)
+{
+    x = rhs.x;
+    y = rhs.y;
+    z = rhs.z;
+    w = rhs.w;
+    return;
+}
+
+// vec4f assignment and comparison operators
+vec4f& vec4f::operator=(const vec4f &rhs)
+{
+    this->x = rhs.x;
+    this->y = rhs.y;
+    this->z = rhs.z;
+    this->w = rhs.w;
+    return *this;
+}
+
+bool vec4f::operator==(vec4f rhs)
+{
+    return ((this->x == rhs.x) && (this->y == rhs.y) && (this->z == rhs.z) && (this->w == rhs.w));
+}
+
+// vec4f arithmetic operators
+vec4f vec4f::operator+(vec4f rhs) const
+{
+    vec4f ret;
+    ret.x = this->x + rhs.x;
+    ret.y = this->y + rhs.y;
+    ret.z = this->z + rhs.z;
+    ret.w = this->w + rhs.w;
+    return ret;
+}
+
+vec4f vec4f::operator-(vec4f rhs) const
+{
+    vec4f ret;
+    ret.x = this->x - rhs.x;
+    ret.y = this->y - rhs.y;
+    ret.z = this->z - rhs.z;
+    ret.w = this->w - rhs.w;
+    return ret;
+}
+
+vec4f vec4f::operator*(float rhs) const
+{
+    vec4f ret;
+    ret.x = this->x*rhs;
+    ret.y = this->y*rhs;
+    ret.z = this->z*rhs;
+    ret.w = this->w*rhs;
+    return ret;
+}
+
+vec4f vec4f::operator/(float rhs) const
+{
+    vec4f ret;
+    ret.x = this->x/rhs;
+    ret.y = this->y/rhs;
+    ret.z = this->z/rhs;
+    ret.w = this->w/rhs;
+    return ret;
+}
+
+vec4f operator-(vec4f rhs)
+{
+    vec4f ret;
+    ret.x = -rhs.x;
+    ret.y = -rhs.y;
+    ret.z = -rhs.z;
+    ret.w = -rhs.w;
+    return ret;
+}
+
+// vec4f subscript operators
+float& vec4f::operator[](int ind)
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    else if(ind==2) return z;
+    else if(ind==3) return w;
+    return x;
+}
+
+const float& vec4f::operator[](int ind) const
+{
+    if(ind==0) return x;
+    else if(ind==1) return y;
+    else if(ind==2) return z;
+    else if(ind==3) return w;
+    return x;
+}
+
+// vector operations (new nomenclature)
+float vec4f_dot(vec4f v1, vec4f v2)
+{
+    float ret = 0.0f;
+    ret += v1.x*v2.x;
+    ret += v1.y*v2.y;
+    ret += v1.z*v2.z;
+    ret += v1.w*v2.w;
+    return ret;
+}
+
+float vec4f_norm(vec4f v)
+{
+    float ret = ::sqrtf(v.x*v.x+v.y*v.y+v.z*v.z+v.w*v.w);
+    return ret;
+}
+
+vec4f vec4f_normalize(vec4f v)
+{
+    vec4f ret = vec4f(0,0,0,0);
+    float nrm = vec4f_norm(v);
+    if(nrm>0.0)
+    {
+        ret = v/nrm;
+    }
+    return ret;
+}
+
+float vec4f_dist(vec4f v1, vec4f v2)
+{
+    return vec4f_norm(v1-v2);
+}
+
+vec4f vec4f_proj(vec4f v1, vec4f v2)
+{
+    vec4f ret;
+    vec4f v2n = vec4f_normalize(v2);
+    ret = v2n*vec4f_dot(v1,v2n);
+    return ret;
+}
+
+// vector operations (old nomenclature)
+float vdot(vec4f v1, vec4f v2)
+{
+    return vec4f_dot(v1,v2);
+}
+
+float vnorm(vec4f v)
+{
+    return vec4f_norm(v);
+}
+
+vec4f vnormalize(vec4f v)
+{
+    return vec4f_normalize(v);
+}
+
+float vdist(vec4f v1, vec4f v2)
+{
+    return vec4f_dist(v1,v2);
+}
+
+vec4f vproj(vec4f v1, vec4f v2)
+{
+    return vec4f_proj(v1,v2);
+}
+
+// mat4f constructors
+mat4f::mat4f()
+{
+    int I;
+    for(I=0;I<16;I++) data[I] = 0.0f;
+    return;
+}
+
+mat4f::mat4f(
+    float _xx, float _yx, float _zx, float _wx,
+    float _xy, float _yy, float _zy, float _wy,
+    float _xz, float _yz, float _zz, float _wz,
+    float _xw, float _yw, float _zw, float _ww
+)
+{
+    data[0] = _xx;
+    data[1] = _yx;
+    data[2] = _zx;
+    data[3] = _wx;
+    
+    data[4] = _xy;
+    data[5] = _yy;
+    data[6] = _zy;
+    data[7] = _wy;
+    
+    data[8] = _xz;
+    data[9] = _yz;
+    data[10] = _zz;
+    data[11] = _wz;
+    
+    data[12] = _xw;
+    data[13] = _yw;
+    data[14] = _zw;
+    data[15] = _ww;
+    
+    return;
+}
+
+mat4f::mat4f(const float *_data)
+{
+    int I;
+    for(I=0;I<16;I++) data[I] = _data[I];
+    return;
+}
+
+mat4f::mat4f(const mat4f &rhs)
+{
+    int I;
+    for(I=0;I<16;I++) data[I] = rhs.data[I];
+    return;
+}
+
+mat4f::mat4f(const vec4f _col1, const vec4f _col2, const vec4f _col3, const vec4f _col4)
+{
+    data[0] = _col1[0];
+    data[1] = _col1[1];
+    data[2] = _col1[2];
+    data[3] = _col1[3];
+
+    data[4] = _col2[0];
+    data[5] = _col2[1];
+    data[6] = _col2[2];
+    data[7] = _col2[3];
+
+    data[8] = _col3[0];
+    data[9] = _col3[1];
+    data[10] = _col3[2];
+    data[11] = _col3[3];
+
+    data[12] = _col4[0];
+    data[13] = _col4[1];
+    data[14] = _col4[2];
+    data[15] = _col4[3];
+    
+    return;
+}
+
+// mat4f assignment and comparison operators
+mat4f& mat4f::operator=(const mat4f rhs)
+{
+    int I;
+    for(I=0;I<16;I++) data[I] = rhs.data[I];
+    return *this;
+}
+
+bool mat4f::operator==(const mat4f rhs)
+{
+    bool ret = 1;
+    int I;
+    for(I=0;I<16;I++)
+    {
+        if(this->data[I]!=rhs.data[I])
+        {
+            ret = 0;
+            break;
+        }
+    }
+
+    return ret;
+}
+
+// mat4f arithmetic operators
+mat4f mat4f::operator+(mat4f rhs) const
+{
+    mat4f ret;
+    int I;
+    for(I=0;I<16;I++)
+    {
+        ret.data[I] = this->data[I] + rhs.data[I];
+    }
+    return ret;
+}
+
+mat4f mat4f::operator-(mat4f rhs) const
+{
+    mat4f ret;
+    int I;
+    for(I=0;I<16;I++)
+    {
+        ret.data[I] = this->data[I] - rhs.data[I];
+    }
+    return ret;
+}
+
+mat4f mat4f::operator*(mat4f rhs) const
+{
+    mat4f ret = mat4f_zeros();
+    int I,J,K;
+    for(I=0;I<4;I++)
+    for(J=0;J<4;J++)
+    for(K=0;K<4;K++)
+        ret(I,K) += (*this)(I,J)*rhs(J,K);
+    return ret;
+}
+
+mat4f mat4f::operator/(mat4f rhs) const
+{
+    mat4f rhsi = rhs.inverse();
+    return (*this)*rhsi;
+}
+
+mat4f mat4f::operator*(float rhs) const
+{
+    mat4f ret;
+    int I;
+    for(I=0;I<16;I++) ret.data[I] = this->data[I]*rhs;
+    return ret;
+}
+
+mat4f mat4f::operator/(float rhs) const
+{
+    mat4f ret;
+    int I;
+    for(I=0;I<16;I++) ret.data[I] = this->data[I]/rhs;
+    return ret;
+}
+
+vec4f mat4f::operator*(vec4f rhs) const
+{
+    vec4f ret = vec4f(0,0,0,0);
+    int I,J;
+    for(I=0;I<4;I++)
+    for(J=0;J<4;J++)
+        ret[I] += (*this)(I,J)*rhs[J];
+    return ret;
+}
+
+// mat4f access operators
+float& mat4f::operator()(int I, int J)
+{
+    int ind = I + 4*J;
+    return data[ind];
+}
+
+float& mat4f::operator[](int I)
+{
+    return data[I];
+}
+
+const float& mat4f::operator()(int I, int J) const
+{
+    int ind = I + 4*J;
+    return data[ind];
+}
+
+const float& mat4f::operator[](int I) const
+{
+    return data[I];
+}
+
+float& mat4f::at(int I, int J)
+{
+    int ind = I + 4*J;
+    return data[ind];
+}
+
+const float& mat4f::at(int I, int J) const
+{
+    int ind = I + 4*J;
+    return data[ind];
+}
+
+float mat4f::det()
+{
+    float ret = 0.0f;
+    float a00,a01,a02,a03;
+    float a10,a11,a12,a13;
+    float a20,a21,a22,a23;
+    float a30,a31,a32,a33;
+  
+    a00 = this->at(0,0);
+    a01 = this->at(0,1);
+    a02 = this->at(0,2);
+    a03 = this->at(0,3);
+    a10 = this->at(1,0);
+    a11 = this->at(1,1);
+    a12 = this->at(1,2);
+    a13 = this->at(1,3);
+    a20 = this->at(2,0);
+    a21 = this->at(2,1);
+    a22 = this->at(2,2);
+    a23 = this->at(2,3);
+    a30 = this->at(3,0);
+    a31 = this->at(3,1);
+    a32 = this->at(3,2);
+    a33 = this->at(3,3);
+  
+    ret = a03*a12*a21*a30 -
+          a02*a13*a21*a30 -
+          a03*a11*a22*a30 +
+          a01*a13*a22*a30 +
+          a02*a11*a23*a30 -
+          a01*a12*a23*a30 -
+          a03*a12*a20*a31 +
+          a02*a13*a20*a31 +
+          a03*a10*a22*a31 -
+          a00*a13*a22*a31 -
+          a02*a10*a23*a31 +
+          a00*a12*a23*a31 +
+          a03*a11*a20*a32 -
+          a01*a13*a20*a32 -
+          a03*a10*a21*a32 +
+          a00*a13*a21*a32 +
+          a01*a10*a23*a32 -
+          a00*a11*a23*a32 -
+          a02*a11*a20*a33 +
+          a01*a12*a20*a33 +
+          a02*a10*a21*a33 -
+          a00*a12*a21*a33 -
+          a01*a10*a22*a33 +
+          a00*a11*a22*a33;
+
+    return ret;
+}
+
+mat4f mat4f::inverse()
+{
+    mat4f ret = mat4f_zeros();
+    float dt = this->det();
+    
+    float a00,a01,a02,a03;
+    float a10,a11,a12,a13;
+    float a20,a21,a22,a23;
+    float a30,a31,a32,a33;
+  
+    float b00,b01,b02,b03;
+    float b10,b11,b12,b13;
+    float b20,b21,b22,b23;
+    float b30,b31,b32,b33;
+  
+    a00 = this->at(0,0);
+    a01 = this->at(0,1);
+    a02 = this->at(0,2);
+    a03 = this->at(0,3);
+    a10 = this->at(1,0);
+    a11 = this->at(1,1);
+    a12 = this->at(1,2);
+    a13 = this->at(1,3);
+    a20 = this->at(2,0);
+    a21 = this->at(2,1);
+    a22 = this->at(2,2);
+    a23 = this->at(2,3);
+    a30 = this->at(3,0);
+    a31 = this->at(3,1);
+    a32 = this->at(3,2);
+    a33 = this->at(3,3);
+  
+    if(ams::abs(dt)>vecmat_det_small)
+    {
+      b00 = -a13*a22*a31 + a12*a23*a31 + a13*a21*a32 - a11*a23*a32 - a12*a21*a33 + a11*a22*a33;
+      b01 =  a03*a22*a31 - a02*a23*a31 - a03*a21*a32 + a01*a23*a32 + a02*a21*a33 - a01*a22*a33;
+      b02 = -a03*a12*a31 + a02*a13*a31 + a03*a11*a32 - a01*a13*a32 - a02*a11*a33 + a01*a12*a33;
+      b03 =  a03*a12*a21 - a02*a13*a21 - a03*a11*a22 + a01*a13*a22 + a02*a11*a23 - a01*a12*a23;
+      b10 =  a13*a22*a30 - a12*a23*a30 - a13*a20*a32 + a10*a23*a32 + a12*a20*a33 - a10*a22*a33;
+      b11 = -a03*a22*a30 + a02*a23*a30 + a03*a20*a32 - a00*a23*a32 - a02*a20*a33 + a00*a22*a33;
+      b12 =  a03*a12*a30 - a02*a13*a30 - a03*a10*a32 + a00*a13*a32 + a02*a10*a33 - a00*a12*a33;
+      b13 = -a03*a12*a20 + a02*a13*a20 + a03*a10*a22 - a00*a13*a22 - a02*a10*a23 + a00*a12*a23;
+      b20 = -a13*a21*a30 + a11*a23*a30 + a13*a20*a31 - a10*a23*a31 - a11*a20*a33 + a10*a21*a33;
+      b21 =  a03*a21*a30 - a01*a23*a30 - a03*a20*a31 + a00*a23*a31 + a01*a20*a33 - a00*a21*a33;
+      b22 = -a03*a11*a30 + a01*a13*a30 + a03*a10*a31 - a00*a13*a31 - a01*a10*a33 + a00*a11*a33;
+      b23 =  a03*a11*a20 - a01*a13*a20 - a03*a10*a21 + a00*a13*a21 + a01*a10*a23 - a00*a11*a23;
+      b30 =  a12*a21*a30 - a11*a22*a30 - a12*a20*a31 + a10*a22*a31 + a11*a20*a32 - a10*a21*a32;
+      b31 = -a02*a21*a30 + a01*a22*a30 + a02*a20*a31 - a00*a22*a31 - a01*a20*a32 + a00*a21*a32;
+      b32 =  a02*a11*a30 - a01*a12*a30 - a02*a10*a31 + a00*a12*a31 + a01*a10*a32 - a00*a11*a32;
+      b33 = -a02*a11*a20 + a01*a12*a20 + a02*a10*a21 - a00*a12*a21 - a01*a10*a22 + a00*a11*a22;
+      b00 = b00/dt;
+      b01 = b01/dt;
+      b02 = b02/dt;
+      b03 = b03/dt;
+      b10 = b10/dt;
+      b11 = b11/dt;
+      b12 = b12/dt;
+      b13 = b13/dt;
+      b20 = b20/dt;
+      b21 = b21/dt;
+      b22 = b22/dt;
+      b23 = b23/dt;
+      b30 = b30/dt;
+      b31 = b31/dt;
+      b32 = b32/dt;
+      b33 = b33/dt;
+      ret(0,0) = b00;
+      ret(0,1) = b01;
+      ret(0,2) = b02;
+      ret(0,3) = b03;
+      ret(1,0) = b10;
+      ret(1,1) = b11;
+      ret(1,2) = b12;
+      ret(1,3) = b13;
+      ret(2,0) = b20;
+      ret(2,1) = b21;
+      ret(2,2) = b22;
+      ret(2,3) = b23;
+      ret(3,0) = b30;
+      ret(3,1) = b31;
+      ret(3,2) = b32;
+      ret(3,3) = b33;
+    }
+    //this was STUPID. Gaah. Computer algebra system saves the day? I'd be surprised if this didn't end up *slower* than gaussian elimination. Don't do this again!
+
+    return ret;
+}
+
+void mat4f::invert()
+{
+    mat4f m = this->inverse();
+    *this = m;
+}
+
+mat4f mat4f::transpose()
+{
+    mat4f ret;
+    int I,J;
+    for(I=0;I<4;I++)
+    for(J=0;J<4;J++)
+        ret(I,J) = this->at(J,I);
+    return ret;
+}
+
+void mat4f::dotranspose()
+{
+    mat4f m = this->transpose();
+    *this = m;
+}
+
+// matrix operations (standalone functions)
+mat4f mat4f_eye()
+{
+    mat4f ret = mat4f_zeros();
+    ret(0,0) = 1.0f;
+    ret(1,1) = 1.0f;
+    ret(2,2) = 1.0f;
+    ret(3,3) = 1.0f;
+    return ret;
+}
+
+mat4f mat4f_zeros()
+{
+    mat4f ret;
+    // for speed - the constructor already sets data[I] to 0.0
+    // int I;
+    // for(I=0;I<16;I++) ret.data[I] = 0.0;
+    return ret;
+}
+
+mat4f mat4f_ones()
+{
+    mat4f ret;
+    int I;
+    for(I=0;I<16;I++) ret.data[I] = 1.0f;
+    return ret;
+}
+
+mat4f mat4f_mult(mat4f a, mat4f b)
+{
+    return a*b;
+}
+
+vec4f mat4f_mult(vec4f a, mat4f b)
+{
+    vec4f ret = vec4f(0,0,0,0);
+    int I,J;
+    for(I=0;I<4;I++)
+    for(J=0;J<4;J++)
+        ret[J] += a[I] * b(I,J);
+
+    return ret;
+}
+
+vec4f mat4f_mult(mat4f a, vec4f b)
+{
+    return a*b;
+}
+
+float mat4f_det(mat4f a)
+{
+    return a.det();
+}
+
+mat4f mat4f_inverse(mat4f a)
+{
+    return a.inverse();
+}
+
+mat4f mat4f_transpose(mat4f a)
+{
+    return a.transpose();
+}
+
+};
\ No newline at end of file
diff --git a/src/amsmathutil25/math/amsmathutil2t_complex128.cpp b/src/amsmathutil25/math/amsmathutil2t_complex128.cpp
new file mode 100644
index 0000000..602fa18
--- /dev/null
+++ b/src/amsmathutil25/math/amsmathutil2t_complex128.cpp
@@ -0,0 +1,414 @@
+
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+namespace cmp
+{
+
+complex::complex()
+{
+  //initcode
+  r = 0.0;
+  i = 0.0;
+  return;
+}
+
+complex::complex(double _r, double _i)
+{
+  r = _r; i = _i;
+  return;
+}
+
+complex::complex(double _r)
+{
+  r = _r; i = 0.0;
+  return;
+}
+
+complex::complex(const complex &other)
+{
+  r = other.r;
+  i = other.i;
+  return;
+}
+
+complex& complex::operator=(const complex& other)
+{
+  r = other.r;
+  i = other.i;
+  return *this;
+}
+
+complex& complex::operator=(const double other)
+{
+  r = other;
+  i = 0.0;
+  return *this;
+}
+
+double& complex::operator[](int& ind)
+{
+  if(ind==0)
+  {
+    return this->r;
+  }
+  else if(ind==1)
+  {
+    return this->i;
+  }
+  else
+  {
+    printf("complex::operator[] index %d out of bounds\n",ind);
+    exit(0);
+  }
+}
+
+const double& complex::operator[](const int& ind) const
+{
+  if(ind==0)
+  {
+    return this->r;
+  }
+  else if(ind==1)
+  {
+    return this->i;
+  }
+  else
+  {
+    printf("complex::operator[] index %d out of bounds\n",ind);
+    exit(0);
+  }
+}
+
+complex complex::operator+(const complex &z)
+{
+  complex cp;
+  cp.r = r+z.r;
+  cp.i = i+z.i;
+  return cp;
+}
+
+ complex  complex::operator-(const complex &z)
+{
+  complex cp;
+  cp.r = r-z.r;
+  cp.i = i-z.i;
+  return cp;
+}
+
+ complex complex::operator*(const complex &z)
+{
+  complex cp;
+	cp.r = r * z.r - i * z.i;
+	cp.i = r * z.i + i * z.r;
+	return cp;
+}
+
+ complex complex::operator/(const complex &z)
+{
+	complex cp;
+	cp.r = (r*z.r+i*z.i)/(z.r*z.r+z.i*z.i);
+	cp.i = (i*z.r-r*z.i)/(z.r*z.r+z.i*z.i);
+	return cp;
+}
+
+ complex complex::operator+(const double &z)
+{
+	complex cp;
+	cp.r = r + z;
+	cp.i = i + 0.0;
+	return cp;
+}
+
+ complex complex::operator-(const double &z)
+{
+	complex cp;
+	cp.r = r - z;
+	cp.i = i - 0.0;
+	return cp;
+}
+
+ complex complex::operator*(const double &z)
+{
+	complex cp;
+	cp.r = r * z;
+	cp.i = i * z;
+	return cp;
+}
+
+ complex complex::operator/(const double &z)
+{
+	complex cp;
+	cp.r = r/z;
+	cp.i = i/z;
+	return cp;
+}
+
+complex operator-(const complex &z)
+{
+	complex ret;
+	ret.r = -z.r;
+	ret.i = -z.i;
+	return ret;
+}
+
+double complex::arg()
+{
+	double ret = ams::arg(this->r,this->i);
+	return ret;
+}
+
+double complex::mag()
+{
+	return ::sqrt(r*r+i*i);
+}
+
+bool complex::iszero()
+{
+	bool ret = 0;
+	if(r==0.0 && i == 0.0)
+  {
+    ret = 1;
+  }
+	return ret;
+}
+
+bool complex::isreal()
+{
+  bool ret = 0;
+	// bool ret = 0;
+	// if(re>DSMALL||re<-DSMALL)
+	// {
+	// 	ret = 0;
+	// }
+	// else
+	// {
+	// 	ret = 1;
+	// }
+  if(i==0.0)
+  {
+    ret = 1;
+  }
+
+	return ret;
+}
+
+bool complex::isimag()
+{
+  bool ret = 0;
+  if(r==0.0)
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+const complex complex::conj()
+{
+  complex cp;
+  cp.r = r;
+  cp.i = -i;
+  return cp;
+}
+
+double real(complex z)
+{
+  return z.r;
+}
+
+double imag(complex z)
+{
+  return z.i;
+}
+
+complex conj(complex z)
+{
+  complex cp;
+  cp.r = z.r;
+  cp.i = -z.i;
+  return cp;
+}
+
+complex exp(complex z)
+{
+  complex ret;
+  ret.r = ::exp(z.r)*::cos(z.i);
+  ret.i = ::exp(z.r)*::sin(z.i);
+  return ret;
+}
+
+
+complex log(complex z)
+{
+  complex ret;
+  ret.r = ::log(::sqrt(z.r*z.r+z.i*z.i));
+  ret.i = ams::arg(z.r,z.i);
+  return ret;
+}
+
+complex cos(complex z)
+{
+  complex ret = complex(0.0,0.0);
+  complex div = complex(2.0,0.0);
+  complex im1 = complex(0.0,1.0);
+  ret = (exp(im1*z)+exp(im1*(-z)))/div;
+  return ret;
+}
+
+complex sin(complex z)
+{
+  complex ret = complex(0.0,0.0);
+  complex div = complex(0.0,2.0);
+  complex im1 = complex(0.0,1.0);
+  ret = (exp(im1*z)-exp(im1*(-z)))/div;
+  return ret;
+}
+
+complex cosh(complex z)
+{
+  complex ret = complex(0.0,0.0);
+  complex div = complex(2.0,0.0);
+
+  ret = (exp(z)+exp((-z)))/div;
+  return ret;
+}
+
+complex sinh(complex z)
+{
+  complex ret = complex(0.0,0.0);
+  complex div = complex(2.0,0.0);
+
+  ret = (exp(z)-exp(-z))/div;
+  return ret;
+}
+
+complex tan(complex z)
+{
+	return sin(z)/cos(z);
+}
+
+complex tanh(complex z)
+{
+	return sinh(z)/cosh(z);
+}
+
+complex pow(complex z1, complex z2)
+{
+	complex ret;
+	if(z1.mag()>DSMALL)
+		ret= ams::cmp::exp(z2*ams::cmp::log(z1));
+	else
+		ret = 0.0;
+	return ret;
+}
+
+double abs(complex z)
+{
+  return z.mag();
+}
+
+double arg(ams::cmp::complex z)
+{
+  return ams::arg(z.r,z.i);
+}
+
+//comparison operators
+bool complex::operator==(const complex &z)
+{
+  bool ret = 0;
+  if((z.r == this->r) && (z.i == this->i))
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+bool complex::operator!=(const complex &z)
+{
+  bool ret = !(*this == z);
+  return ret;
+}
+
+//sort first by real value, then by imaginary value
+//this is done so that an ordering exists, as long as two values aren't equal
+bool complex::operator>(const complex& z)
+{
+  bool ret = 0;
+  if(this->r > z.r)
+  {
+    ret = 1;
+  }
+  else if(this->r == z.r)
+  {
+    if(this->i > z.i)
+    {
+      ret = 1;
+    }
+  }
+
+  return ret;
+}
+
+bool complex::operator<(const complex& z)
+{
+  bool ret = 0;
+  ret = !((*this>z) || *this==z);
+  return ret;
+}
+
+bool complex::operator>=(const complex &z)
+{
+  bool ret = 0;
+  ret = (*this==z) || (*this>z);
+  return ret;
+}
+
+bool complex::operator<=(const complex &z)
+{
+  bool ret = 0;
+  ret = (*this<z) || (*this==z);
+  return ret;
+}
+
+bool complex::isnan()
+{
+  bool ret = 0;
+  if(ams::isnan(this->r) || ams::isnan(this->i))
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+bool complex::isinf()
+{
+  bool ret = 0;
+  //calls math.h isinf()
+  if(::isinf(this->r) || ::isinf(this->i))
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+complex csgn(complex z)
+{
+  double a;
+  complex ret = ams::cmp::complex(0,0);
+  if(z.r==0.0 && z.i==0.0)
+  {
+    //do nothing
+  }
+  else
+  {
+    a = ams::cmp::arg(z);
+    ret.r = ::cos(a);
+    ret.i = ::sin(a);
+  }
+
+  return ret;
+}
+
+};//end namespace cmp
+};//end namespace ams
\ No newline at end of file
diff --git a/src/amsmathutil25/math/amsmathutil2t_complex64.cpp b/src/amsmathutil25/math/amsmathutil2t_complex64.cpp
new file mode 100644
index 0000000..16fce24
--- /dev/null
+++ b/src/amsmathutil25/math/amsmathutil2t_complex64.cpp
@@ -0,0 +1,414 @@
+
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+namespace cmp
+{
+
+complex64::complex64()
+{
+  //initcode
+  r = 0.0;
+  i = 0.0;
+  return;
+}
+
+complex64::complex64(float _r, float _i)
+{
+  r = _r; i = _i;
+  return;
+}
+
+complex64::complex64(float _r)
+{
+  r = _r; i = 0.0;
+  return;
+}
+
+complex64::complex64(const complex64 &other)
+{
+  r = other.r;
+  i = other.i;
+  return;
+}
+
+complex64& complex64::operator=(const complex64& other)
+{
+  r = other.r;
+  i = other.i;
+  return *this;
+}
+
+complex64& complex64::operator=(const float other)
+{
+  r = other;
+  i = 0.0;
+  return *this;
+}
+
+float& complex64::operator[](int& ind)
+{
+  if(ind==0)
+  {
+    return this->r;
+  }
+  else if(ind==1)
+  {
+    return this->i;
+  }
+  else
+  {
+    printf("complex64::operator[] index %d out of bounds\n",ind);
+    exit(0);
+  }
+}
+
+const float& complex64::operator[](const int& ind) const
+{
+  if(ind==0)
+  {
+    return this->r;
+  }
+  else if(ind==1)
+  {
+    return this->i;
+  }
+  else
+  {
+    printf("complex64::operator[] index %d out of bounds\n",ind);
+    exit(0);
+  }
+}
+
+complex64 complex64::operator+(const complex64 &z)
+{
+  complex64 cp;
+  cp.r = r+z.r;
+  cp.i = i+z.i;
+  return cp;
+}
+
+ complex64  complex64::operator-(const complex64 &z)
+{
+  complex64 cp;
+  cp.r = r-z.r;
+  cp.i = i-z.i;
+  return cp;
+}
+
+ complex64 complex64::operator*(const complex64 &z)
+{
+  complex64 cp;
+	cp.r = r * z.r - i * z.i;
+	cp.i = r * z.i + i * z.r;
+	return cp;
+}
+
+ complex64 complex64::operator/(const complex64 &z)
+{
+	complex64 cp;
+	cp.r = (r*z.r+i*z.i)/(z.r*z.r+z.i*z.i);
+	cp.i = (i*z.r-r*z.i)/(z.r*z.r+z.i*z.i);
+	return cp;
+}
+
+ complex64 complex64::operator+(const float &z)
+{
+	complex64 cp;
+	cp.r = r + z;
+	cp.i = i + 0.0;
+	return cp;
+}
+
+ complex64 complex64::operator-(const float &z)
+{
+	complex64 cp;
+	cp.r = r - z;
+	cp.i = i - 0.0;
+	return cp;
+}
+
+ complex64 complex64::operator*(const float &z)
+{
+	complex64 cp;
+	cp.r = r * z;
+	cp.i = i * z;
+	return cp;
+}
+
+ complex64 complex64::operator/(const float &z)
+{
+	complex64 cp;
+	cp.r = r/z;
+	cp.i = i/z;
+	return cp;
+}
+
+complex64 operator-(const complex64 &z)
+{
+	complex64 ret;
+	ret.r = -z.r;
+	ret.i = -z.i;
+	return ret;
+}
+
+float complex64::arg()
+{
+	float ret = ams::arg(this->r,this->i);
+	return ret;
+}
+
+float complex64::mag()
+{
+	return ::sqrtf(r*r+i*i);
+}
+
+bool complex64::iszero()
+{
+	bool ret = 0;
+	if(r==0.0 && i == 0.0)
+  {
+    ret = 1;
+  }
+	return ret;
+}
+
+bool complex64::isreal()
+{
+  bool ret = 0;
+	// bool ret = 0;
+	// if(re>DSMALL||re<-DSMALL)
+	// {
+	// 	ret = 0;
+	// }
+	// else
+	// {
+	// 	ret = 1;
+	// }
+  if(i==0.0)
+  {
+    ret = 1;
+  }
+
+	return ret;
+}
+
+bool complex64::isimag()
+{
+  bool ret = 0;
+  if(r==0.0)
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+const complex64 complex64::conj()
+{
+  complex64 cp;
+  cp.r = r;
+  cp.i = -i;
+  return cp;
+}
+
+float real(complex64 z)
+{
+  return z.r;
+}
+
+float imag(complex64 z)
+{
+  return z.i;
+}
+
+complex64 conj(complex64 z)
+{
+  complex64 cp;
+  cp.r = z.r;
+  cp.i = -z.i;
+  return cp;
+}
+
+complex64 exp(complex64 z)
+{
+  complex64 ret;
+  ret.r = ::expf(z.r)*::cos(z.i);
+  ret.i = ::expf(z.r)*::sin(z.i);
+  return ret;
+}
+
+
+complex64 log(complex64 z)
+{
+  complex64 ret;
+  ret.r = ::log(::sqrtf(z.r*z.r+z.i*z.i));
+  ret.i = ams::arg(z.r,z.i);
+  return ret;
+}
+
+complex64 cos(complex64 z)
+{
+  complex64 ret = complex64(0.0,0.0);
+  complex64 div = complex64(2.0,0.0);
+  complex64 im1 = complex64(0.0,1.0);
+  ret = (exp(im1*z)+exp(im1*(-z)))/div;
+  return ret;
+}
+
+complex64 sin(complex64 z)
+{
+  complex64 ret = complex64(0.0,0.0);
+  complex64 div = complex64(0.0,2.0);
+  complex64 im1 = complex64(0.0,1.0);
+  ret = (exp(im1*z)-exp(im1*(-z)))/div;
+  return ret;
+}
+
+complex64 cosh(complex64 z)
+{
+  complex64 ret = complex64(0.0,0.0);
+  complex64 div = complex64(2.0,0.0);
+
+  ret = (exp(z)+exp((-z)))/div;
+  return ret;
+}
+
+complex64 sinh(complex64 z)
+{
+  complex64 ret = complex64(0.0,0.0);
+  complex64 div = complex64(2.0,0.0);
+
+  ret = (exp(z)-exp(-z))/div;
+  return ret;
+}
+
+complex64 tan(complex64 z)
+{
+	return sin(z)/cos(z);
+}
+
+complex64 tanh(complex64 z)
+{
+	return sinh(z)/cosh(z);
+}
+
+complex64 pow(complex64 z1, complex64 z2)
+{
+	complex64 ret;
+	if(z1.mag()>DSMALL)
+		ret= ams::cmp::exp(z2*ams::cmp::log(z1));
+	else
+		ret = 0.0;
+	return ret;
+}
+
+float abs(complex64 z)
+{
+  return z.mag();
+}
+
+float arg(ams::cmp::complex64 z)
+{
+  return ams::arg(z.r,z.i);
+}
+
+//comparison operators
+bool complex64::operator==(const complex64 &z)
+{
+  bool ret = 0;
+  if((z.r == this->r) && (z.i == this->i))
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+bool complex64::operator!=(const complex64 &z)
+{
+  bool ret = !(*this == z);
+  return ret;
+}
+
+//sort first by real value, then by imaginary value
+//this is done so that an ordering exists, as long as two values aren't equal
+bool complex64::operator>(const complex64& z)
+{
+  bool ret = 0;
+  if(this->r > z.r)
+  {
+    ret = 1;
+  }
+  else if(this->r == z.r)
+  {
+    if(this->i > z.i)
+    {
+      ret = 1;
+    }
+  }
+
+  return ret;
+}
+
+bool complex64::operator<(const complex64& z)
+{
+  bool ret = 0;
+  ret = !((*this>z) || *this==z);
+  return ret;
+}
+
+bool complex64::operator>=(const complex64 &z)
+{
+  bool ret = 0;
+  ret = (*this==z) || (*this>z);
+  return ret;
+}
+
+bool complex64::operator<=(const complex64 &z)
+{
+  bool ret = 0;
+  ret = (*this<z) || (*this==z);
+  return ret;
+}
+
+bool complex64::isnan()
+{
+  bool ret = 0;
+  if(ams::isnan(this->r) || ams::isnan(this->i))
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+bool complex64::isinf()
+{
+  bool ret = 0;
+  //calls math.h isinf()
+  if(::isinf(this->r) || ::isinf(this->i))
+  {
+    ret = 1;
+  }
+  return ret;
+}
+
+complex64 csgn(complex64 z)
+{
+  float a;
+  complex64 ret = ams::cmp::complex64(0,0);
+  if(z.r==0.0 && z.i==0.0)
+  {
+    //do nothing
+  }
+  else
+  {
+    a = ams::cmp::arg(z);
+    ret.r = ::cos(a);
+    ret.i = ::sin(a);
+  }
+
+  return ret;
+}
+
+};//end namespace cmp
+};//end namespace ams
\ No newline at end of file
diff --git a/src/amsmathutil25/random/amsmathutiil25_random.cpp b/src/amsmathutil25/random/amsmathutiil25_random.cpp
new file mode 100644
index 0000000..de8a401
--- /dev/null
+++ b/src/amsmathutil25/random/amsmathutiil25_random.cpp
@@ -0,0 +1,29 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+namespace rand
+{
+
+amsmu_randt1 dpr32_rseed;
+amsmu_randt2 dpr64_rseed;
+
+amsmu_randt1 dpr32_nextseed(amsmu_randt1 seed)
+{
+    amsmu_randt1 sret = seed;
+    sret = ams::mod(sret*dpr32_mult1+dpr32_add1,dpr32_mod);
+    return sret;
+}
+
+amsmu_randt2 dpr64_nextseed(amsmu_randt2 seed)
+{
+    amsmu_randt2 sret = seed;
+    sret = ams::mod(sret*dpr64_mult1+dpr64_add1,dpr64_mod);
+    return sret;
+}
+
+
+
+
+};
+};
\ No newline at end of file
diff --git a/src/amsmathutil25/testing/amsmathtuil25_test1.cpp b/src/amsmathutil25/testing/amsmathtuil25_test1.cpp
index 696873f..ea9a350 100644
--- a/src/amsmathutil25/testing/amsmathtuil25_test1.cpp
+++ b/src/amsmathutil25/testing/amsmathtuil25_test1.cpp
@@ -92,5 +92,7 @@ void test_amsarray2()
 }
 
 
+
+
 }; //end namespace amsmathutil25
 }; //end namespace ams
\ No newline at end of file
diff --git a/src/amsmathutil25/util/amsmathutil25_amsarray.cpp b/src/amsmathutil25/util/amsmathutil25_amsarray.cpp
index 663a136..15a6c06 100644
--- a/src/amsmathutil25/util/amsmathutil25_amsarray.cpp
+++ b/src/amsmathutil25/util/amsmathutil25_amsarray.cpp
@@ -79,6 +79,10 @@ template<> void amsarray<double>::print(bool newline,int printstyle)
     return;
 }
 
-
+// Explicit Class Instantiations?
+template class amsarray<int>; 
+template class amsarray<float>; 
+template class amsarray<int64_t>; 
+template class amsarray<double>; 
 
 };
\ No newline at end of file
diff --git a/src/amsmathutil25/util/amsmathutil25_amsarray_sort.cpp b/src/amsmathutil25/util/amsmathutil25_amsarray_sort.cpp
index 8b7af57..ed5add2 100644
--- a/src/amsmathutil25/util/amsmathutil25_amsarray_sort.cpp
+++ b/src/amsmathutil25/util/amsmathutil25_amsarray_sort.cpp
@@ -59,4 +59,13 @@ amsarray<amsarray_size_t> permutation_identity(amsarray_size_t _length)
     return ret;
 }
 
+void amsarray_permutation_swap(amsarray<amsarray_size_t> *permutation, amsarray_size_t I, amsarray_size_t J)
+{
+    amsarray_size_t tmp;
+    tmp = permutation->data[I];
+    permutation->data[J] = permutation->data[I];
+    permutation->data[I] = tmp;
+    return;
+}
+
 };
\ No newline at end of file