updates

include/amsmathutil25/amsmathutil25.hpp

@@ -9,6 +9,7 @@
 #include <vector>
 #include <thread>
 #include <functional>
+#include <mutex>
 
 namespace ams
 {

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

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

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
 

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
 

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

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
 

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

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

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
 

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
 

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
 

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
 

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();
 

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
+        {
+            for(I=0;I<length && I<other.size();I++)
             {
-                return *this;
-            }
-            else
-            {
-                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
+}
 
 
 

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;
+}
 
 
 

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>

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