diff --git a/amsmathutil25.code-workspace b/amsmathutil25.code-workspace
index 8f4183e..a39253d 100644
--- a/amsmathutil25.code-workspace
+++ b/amsmathutil25.code-workspace
@@ -2,6 +2,12 @@
 	"folders": [
 		{
 			"path": "."
+		},
+		{
+			"path": "../../sourceprojs23/amsmathutil2"
+		},
+		{
+			"path": "../amscppnarray"
 		}
 	]
 }
diff --git a/build/make.linux64.lib.py b/build/make.linux64.lib.py
index 8014d2a..1f915db 100644
--- a/build/make.linux64.lib.py
+++ b/build/make.linux64.lib.py
@@ -16,7 +16,7 @@ builddir = "./build_linux64"
 
 doinstall = False   #copies the build_output to the install dir when finished
 cc = "g++" #compiler
-cflags = "-fPIC"
+cflags = "-O3 -fPIC"
 libraries = "-l{}".format(libname)
 libdirs = "-L{} -L{}/lib -L{}/lib".format(builddir,commondir,depdir)
 linkerflags = "-static -static-libgcc -Wl,-rpath=."
diff --git a/build/make.linux64.test.py b/build/make.linux64.test.py
index d8c1b14..7fef03e 100644
--- a/build/make.linux64.test.py
+++ b/build/make.linux64.test.py
@@ -16,7 +16,7 @@ builddir = "./build_linux64"
 
 doinstall = False   #copies the build_output to the install dir when finished
 cc = "g++" #compiler
-cflags = "-fPIC"
+cflags = "-O3 -fPIC"
 libraries = "-l{}".format(libname)
 libdirs = "-L{} -L{}/lib -L{}/lib".format(builddir,commondir,depdir)
 linkerflags = "-static -static-libgcc -Wl,-rpath=."
diff --git a/build_linux64/libamsmathutil25.linux64.a b/build_linux64/libamsmathutil25.linux64.a
index c7a6dc5..30a862c 100644
Binary files a/build_linux64/libamsmathutil25.linux64.a and b/build_linux64/libamsmathutil25.linux64.a differ
diff --git a/build_linux64/objstore/amsmathtuil25_test1.o b/build_linux64/objstore/amsmathtuil25_test1.o
new file mode 100644
index 0000000..2fbeb00
Binary files /dev/null and b/build_linux64/objstore/amsmathtuil25_test1.o differ
diff --git a/build_linux64/objstore/amsmathutil25_amsarray.o b/build_linux64/objstore/amsmathutil25_amsarray.o
new file mode 100644
index 0000000..aefbb24
Binary files /dev/null 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
new file mode 100644
index 0000000..389530e
Binary files /dev/null and b/build_linux64/objstore/amsmathutil25_amsarray_sort.o differ
diff --git a/build_linux64/tests b/build_linux64/tests
index 125e77e..1b1d8d6 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 5f927ab..c059954 100644
--- a/include/amsmathutil25/amsmathutil25.hpp
+++ b/include/amsmathutil25/amsmathutil25.hpp
@@ -4,11 +4,37 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include <stdint.h>
 
+#include <vector>
+#include <thread>
+#include <functional>
 
 namespace ams
 {
 
+//return values used for this library
+static const int amsmathutil25_success = 1;
+static const int amsmathutil25_failure = -1;
+
+//maximum number of threads to use
+static const int amsmathutil25_maxthreads = 50;
+
+//problem size at which to begin using multi-threading
+static const int amsmathutil25_threadpsz = 5000;
+
+
 };
 
+//Library subsections
+#include <amsmathutil25/math/amsmathutil25_math.hpp>
+#include <amsmathutil25/util/amsmathutil25_util.hpp>
+#include <amsmathutil25/random/amsmathutil25_random.hpp>
+
+
+//... 
+// Testing
+#include <amsmathutil25/testing/amsmathutil25_testing.hpp>
+
+
 #endif
\ No newline at end of file
diff --git a/include/amsmathutil25/math/amsmathutil25_complex128.hpp b/include/amsmathutil25/math/amsmathutil25_complex128.hpp
new file mode 100644
index 0000000..54e9c78
--- /dev/null
+++ b/include/amsmathutil25/math/amsmathutil25_complex128.hpp
@@ -0,0 +1,14 @@
+#ifndef __AMSMATHUTIL25_COMPLEX128_HPP__
+#define __AMSMATHUTIL25_COMPLEX128_HPP__
+
+namespace ams
+{
+namespace cmp
+{
+
+    
+}; //end namespace cmp
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/math/amsmathutil25_complex64.hpp b/include/amsmathutil25/math/amsmathutil25_complex64.hpp
new file mode 100644
index 0000000..cd1f22a
--- /dev/null
+++ b/include/amsmathutil25/math/amsmathutil25_complex64.hpp
@@ -0,0 +1,14 @@
+#ifndef __AMSMATHUTIL25_COMPLEX64_HPP__
+#define __AMSMATHUTIL25_COMPLEX64_HPP__
+
+namespace ams
+{
+namespace cmp
+{
+
+    
+}; //end namespace cmp
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/math/amsmathutil25_math.hpp b/include/amsmathutil25/math/amsmathutil25_math.hpp
new file mode 100644
index 0000000..c64d798
--- /dev/null
+++ b/include/amsmathutil25/math/amsmathutil25_math.hpp
@@ -0,0 +1,22 @@
+#ifndef __AMSMATHUTIL25_MATH_HPP__
+#define __AMSMATHUTIL25_MATH_HPP__
+
+namespace ams
+{
+
+
+}; //end namespace ams
+
+#include <amsmathutil25/math/amsmathutil25_mathfns1.hpp>
+#include <amsmathutil25/math/amsmathutil25_complex64.hpp>
+#include <amsmathutil25/math/amsmathutil25_complex128.hpp>
+#include <amsmathutil25/math/amsmathutil25_vec2.hpp>
+#include <amsmathutil25/math/amsmathutil25_vec3.hpp>
+#include <amsmathutil25/math/amsmathutil25_vec4.hpp>
+#include <amsmathutil25/math/amsmathutil25_vec2f.hpp>
+#include <amsmathutil25/math/amsmathutil25_vec3f.hpp>
+#include <amsmathutil25/math/amsmathutil25_vec4f.hpp>
+
+
+#endif
+
diff --git a/include/amsmathutil25/math/amsmathutil25_mathfns1.hpp b/include/amsmathutil25/math/amsmathutil25_mathfns1.hpp
new file mode 100644
index 0000000..27c6cb4
--- /dev/null
+++ b/include/amsmathutil25/math/amsmathutil25_mathfns1.hpp
@@ -0,0 +1,11 @@
+#ifndef __AMSMATHUTIL25_MATHFNS1_HPP__
+#define __AMSMATHUTIL25_MATHFNS1_HPP__
+
+namespace ams
+{
+
+
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/math/amsmathutil25_vec2.hpp b/include/amsmathutil25/math/amsmathutil25_vec2.hpp
new file mode 100644
index 0000000..189b42f
--- /dev/null
+++ b/include/amsmathutil25/math/amsmathutil25_vec2.hpp
@@ -0,0 +1,36 @@
+#ifndef __AMSMATHUTIL25_VEC2_HPP__
+#define __AMSMATHUTIL25_VEC2_HPP__
+
+namespace ams
+{
+
+class vec2
+{
+public:
+    double x;
+    double y;
+
+    vec2();
+    vec2(double _x, double _y);
+    vec2 operator=(vec2 rhs);
+
+
+};
+
+class mat2
+{
+public:
+    double data[4];
+
+    mat2();
+    mat2(double _xx, double _xy, double _yx, double _yy);
+    mat2(const double *_data);
+    mat2(mat2 &rhs);
+
+    
+};
+
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/math/amsmathutil25_vec2f.hpp b/include/amsmathutil25/math/amsmathutil25_vec2f.hpp
new file mode 100644
index 0000000..00ef248
--- /dev/null
+++ b/include/amsmathutil25/math/amsmathutil25_vec2f.hpp
@@ -0,0 +1,11 @@
+#ifndef __AMSMATHUTIL25_VEC2F_HPP__
+#define __AMSMATHUTIL25_VEC2F_HPP__
+
+namespace ams
+{
+
+
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/math/amsmathutil25_vec3.hpp b/include/amsmathutil25/math/amsmathutil25_vec3.hpp
new file mode 100644
index 0000000..3fca66e
--- /dev/null
+++ b/include/amsmathutil25/math/amsmathutil25_vec3.hpp
@@ -0,0 +1,11 @@
+#ifndef __AMSMATHUTIL25_VEC3_HPP__
+#define __AMSMATHUTIL25_VEC3_HPP__
+
+namespace ams
+{
+
+
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/math/amsmathutil25_vec3f.hpp b/include/amsmathutil25/math/amsmathutil25_vec3f.hpp
new file mode 100644
index 0000000..4d95651
--- /dev/null
+++ b/include/amsmathutil25/math/amsmathutil25_vec3f.hpp
@@ -0,0 +1,11 @@
+#ifndef __AMSMATHUTIL25_VEC3F_HPP__
+#define __AMSMATHUTIL25_VEC3F_HPP__
+
+namespace ams
+{
+
+
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/math/amsmathutil25_vec4.hpp b/include/amsmathutil25/math/amsmathutil25_vec4.hpp
new file mode 100644
index 0000000..21264e1
--- /dev/null
+++ b/include/amsmathutil25/math/amsmathutil25_vec4.hpp
@@ -0,0 +1,11 @@
+#ifndef __AMSMATHUTIL25_VEC4_HPP__
+#define __AMSMATHUTIL25_VEC4_HPP__
+
+namespace ams
+{
+
+
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/math/amsmathutil25_vec4f.hpp b/include/amsmathutil25/math/amsmathutil25_vec4f.hpp
new file mode 100644
index 0000000..631840e
--- /dev/null
+++ b/include/amsmathutil25/math/amsmathutil25_vec4f.hpp
@@ -0,0 +1,11 @@
+#ifndef __AMSMATHUTIL25_VEC4F_HPP__
+#define __AMSMATHUTIL25_VEC4F_HPP__
+
+namespace ams
+{
+
+
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/random/amsmathutil25_random.hpp b/include/amsmathutil25/random/amsmathutil25_random.hpp
new file mode 100644
index 0000000..9389861
--- /dev/null
+++ b/include/amsmathutil25/random/amsmathutil25_random.hpp
@@ -0,0 +1,13 @@
+#ifndef __AMSMATHUTIL25_RANDOM_HPP__
+#define __AMSMATHUTIL25_RANDOM_HPP__
+
+namespace ams
+{
+namespace rand
+{
+
+}; //end namespace rand
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/testing/amsmathutil25_testing.hpp b/include/amsmathutil25/testing/amsmathutil25_testing.hpp
new file mode 100644
index 0000000..7ffaa92
--- /dev/null
+++ b/include/amsmathutil25/testing/amsmathutil25_testing.hpp
@@ -0,0 +1,16 @@
+#ifndef __AMSMATHUTIL25_TESTING_HPP__
+#define __AMSMATHUTIL25_TESTING_HPP__
+
+namespace ams
+{
+namespace amsmathutil25
+{
+
+    void test_amsarray1();
+    void test_amsarray2();
+
+}; //end namespace amsmathutil25
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/util/amsmathutil25_amsarray.hpp b/include/amsmathutil25/util/amsmathutil25_amsarray.hpp
new file mode 100644
index 0000000..c3c89c3
--- /dev/null
+++ b/include/amsmathutil25/util/amsmathutil25_amsarray.hpp
@@ -0,0 +1,116 @@
+#ifndef __AMSMATHUTIL25_AMSARRAY_HPP__
+#define __AMSMATHUTIL25_AMSARRAY_HPP__
+
+namespace ams
+{
+    typedef int64_t amsarray_size_t;
+    static const int amsarray_success = amsmathutil25_success;
+    static const int amsarray_failure = amsmathutil25_failure;
+
+    template<typename T> class amsarray
+    {
+        public:
+        amsarray_size_t length;
+        T *data;
+
+        //Rule of 5 boilerplate
+        amsarray();
+        ~amsarray();
+        amsarray(const amsarray<T>& other);
+        amsarray(amsarray<T>&& other) noexcept;
+        amsarray<T>& operator=(const amsarray<T> &other);
+        amsarray<T>& operator=(amsarray<T> &&other) noexcept;
+
+        //other initializers
+        amsarray<T>& operator=(const std::vector<T>& other);
+        amsarray(const std::initializer_list<T> initlist);
+        amsarray<T>& operator=(const std::initializer_list<T> initlist);
+
+        //casting datatypes
+        template<typename T2> operator amsarray<T2>() const; //casting datatypes
+
+        //resizing operator
+        // returns: 
+        //      amsarray_success on successfully allocating new memory, else
+        //      amsarray_failure.
+        int resize(amsarray_size_t _newlen);
+
+        //returns the array length for routines that expect size() to be present
+        const amsarray_size_t size() const;
+
+        T& operator[](amsarray_size_t ind);
+        const T& operator[](amsarray_size_t ind) const;
+        T& at(amsarray_size_t ind);
+        const T& at(amsarray_size_t ind) const;
+
+        int setall(const T &val);
+
+        //Comparators
+        bool operator==(const amsarray<T>& other) const;
+        bool operator!=(const amsarray<T>& other) const;
+
+        //inserts a value at (before) index ind
+        //insert(0,val) would give {val, oldval0, oldval1, ...}
+        //insert(N,val)
+        int insert(amsarray_size_t ind, const T& val);
+        
+        //removes the value at index ind, and reduces array size by 1
+        int erase(amsarray_size_t ind);
+
+        //aliases to insert/erase operations
+        int append(const T& val);
+        int prepend(const T& val);
+        int push_back(const T& val);
+        int push_front(const T& val);
+        T pop_back();
+        T pop_front();
+
+        //finds the first instance of val in the array
+        amsarray_size_t find(const T& val);
+
+        //finds the next instance of val, starting consideration at indstart
+        amsarray_size_t findnext(amsarray_size_t indstart, const T& val);
+
+        //finds where to insert a particular value such that the list remains ordered
+        amsarray_size_t find_insert_ordered(const T& val);
+
+        //subarray operations
+        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);
+        
+        //returns an array of indices that is a permutation which will sort
+        //this array in ascending order
+        amsarray<amsarray_size_t> sort_permutation();
+
+        //returns an array that is this array in reverse order
+        amsarray<T> reverse();
+
+
+        //C++11 iterators
+        // (magic boilerplate I don't really understand yet)
+        typedef T* iterator;
+        typedef const T* const_iterator;
+        iterator begin() {return data;}
+        iterator end() {return data + length;}
+        const_iterator begin() const {return data;}
+        const_iterator end() const {return data + length;}
+        const_iterator cbegin() const { return data; }
+        const_iterator cend() const { return data+length;}
+
+        //template specialized for particular array types
+        void print(bool newline=0, int printstyle=0);
+    };
+
+// {0, 1, 2, .... N}
+amsarray<amsarray_size_t> permutation_identity(amsarray_size_t _length);
+
+
+}; //end namespace ams
+
+#include <amsmathutil25/util/amsmathutil25_amsarray_impl.hpp>
+#include <amsmathutil25/util/amsmathutil25_amsarray_sortimpl.hpp>
+
+#endif
+
diff --git a/include/amsmathutil25/util/amsmathutil25_amsarray_impl.hpp b/include/amsmathutil25/util/amsmathutil25_amsarray_impl.hpp
new file mode 100644
index 0000000..70fc738
--- /dev/null
+++ b/include/amsmathutil25/util/amsmathutil25_amsarray_impl.hpp
@@ -0,0 +1,747 @@
+#ifndef __AMSMATHUTIL25_AMSARRAY_IMPL_HPP__
+#define __AMSMATHUTIL25_AMSARRAY_IMPL_HPP__
+
+namespace ams
+{
+
+    template<typename T> amsarray<T>::amsarray()
+    {
+        length = 0;
+        data = NULL;
+        return;
+    }
+
+    template<typename T> amsarray<T>::~amsarray()
+    {
+        length = 0;
+        if(data!=NULL) {delete[] data; data = NULL;}
+        return;
+    }
+
+    template<typename T> const amsarray_size_t amsarray<T>::size() const
+    {
+        return this->length;
+    }
+
+    template<typename T> int amsarray<T>::resize(amsarray_size_t _newlen)
+    {
+        int ret = amsarray_success;
+        T *newdata = NULL;
+        amsarray_size_t lmin;
+        T defval = T();
+
+        if(_newlen<=0)
+        {
+            length = 0;
+            if(data!=NULL) {delete[] data; data = NULL;}
+            ret = amsarray_success;
+            return ret;
+        }
+
+        newdata = new(std::nothrow) T[_newlen];
+        if(newdata==NULL)
+        {
+            ret = amsarray_failure;
+            return ret;
+        }
+
+        if(data!=NULL)
+        {
+            lmin = (_newlen>=length) ? length : _newlen;
+            ams::buffer_cast_copy<T,T>(newdata,data,lmin);
+        }
+        ams::buffer_set<T>(newdata,length,_newlen,defval);
+
+        if(data!=NULL) {delete[] data; data = NULL;}
+        data = newdata;
+        length = _newlen;
+
+        return ret;
+    }
+
+    template<typename T> amsarray<T>::amsarray(const amsarray<T>& other)
+    {
+        int res = amsarray_success;
+        length = 0;
+        data = NULL;
+        if(this!=&other)
+        {
+            res = this->resize(other.length);
+            if(res==amsarray_success)
+            {
+                buffer_cast_copy<T,T>(this->data,other.data,length);
+            }
+        }
+        return;
+    }
+
+    template<typename T> amsarray<T>::amsarray(amsarray<T>&& other) noexcept
+    {
+        length = 0;
+        data = NULL;
+        if(this!=&other)
+        {
+            if(data!=NULL) {delete[] data; data = NULL;}
+            length = other.length;
+            data = other.data;
+            other.length = 0;
+            other.data = NULL;
+        }
+        return;
+    }
+
+    template<typename T> amsarray<T>& amsarray<T>::operator=(const amsarray<T> &other)
+    {
+        int res = amsarray_success;
+        if(this!=&other)
+        {
+            res = this->resize(other.length);
+            if(res==amsarray_success)
+            {
+                buffer_cast_copy<T,T>(this->data,other.data,length);
+            }
+        }
+        return *this;
+    }
+
+    template<typename T> amsarray<T>& amsarray<T>::operator=(amsarray<T> &&other) noexcept
+    {
+        if(this!=&other)
+        {
+            if(data!=NULL) {delete[] data; data = NULL;}
+            length = other.length;
+            data = other.data;
+            other.length = 0;
+            other.data = NULL;
+        }
+        return *this;
+    }
+
+    template<typename T> amsarray<T>& amsarray<T>::operator=(const std::vector<T>& other)
+    {
+        amsarray_size_t I;
+        int res;
+    
+        if(this!=&other)
+        {
+            res = this->resize(other.size());
+            if(res!=amsarray_success)
+            {
+                return *this;
+            }
+            else
+            {
+                for(I=0;I<length && I<other.size();I++)
+                {
+                    data[I] = other[I];
+                }
+            }
+        }
+        return *this;
+    }
+    
+    template<typename T> amsarray<T>::amsarray(const std::initializer_list<T> initlist)
+    {
+        amsarray_size_t I;
+        int res;
+        length = 0;
+        data = NULL;
+    
+        res = this->resize(initlist.size());
+        if(res!=amsarray_success)
+        {
+            return;
+        }
+        else
+        {
+            I = 0;
+            for(T elem : initlist)
+            {
+                data[I] = elem;
+                if(I<length) I++;
+            }
+        }
+    }
+    
+    template<typename T> amsarray<T>& amsarray<T>::operator=(const std::initializer_list<T> initlist)
+    {
+        amsarray_size_t I;
+        int res;
+    
+        res = this->resize(initlist.size());
+        if(res!=amsarray_success)
+        {
+            this->length = 0;
+            if(this->data != NULL) {delete[] this->data; this->data=NULL;}
+            return *this;
+        }
+        else
+        {
+            I = 0;
+            for(T elem : initlist)
+            {
+                data[I] = elem;
+                if(I<length) I++;
+            }
+        }
+    
+        return *this;
+    }
+    
+    template<class T> template<class T2> amsarray<T>::operator amsarray<T2>() const
+    {
+        //casting datatypes
+        amsarray<T2> ret;
+        int res;
+    
+        res = ret.resize(this->length);
+        if(res!=amsarray_success)
+        {
+            ret.resize(0);
+            return ret;
+        }
+        else
+        {
+            buffer_cast_copy<T2,T>(ret.data,data,length);
+        }
+    
+        return ret;
+    }
+
+
+
+    template<typename T> T& amsarray<T>::operator[](amsarray_size_t ind)
+    {
+        return data[ind];
+    }
+    template<typename T> const T& amsarray<T>::operator[](amsarray_size_t ind) const
+    {
+        return data[ind];
+    }
+    template<typename T> T& amsarray<T>::at(amsarray_size_t ind)
+    {
+        return data[ind];
+    }
+    template<typename T> const T& amsarray<T>::at(amsarray_size_t ind) const
+    {
+        return data[ind];
+    }
+
+    template<typename T> int amsarray<T>::setall(const T &val)
+    {
+        int ret = amsarray_success;
+        int res;
+        res = buffer_set<T>(data,length,val);
+        if(res!=amsmathutil25_success)
+        {
+            ret = amsarray_failure;
+        }
+        return ret;
+    }
+
+
+
+    template<typename T> amsarray<T> amsarray<T>::subarray(amsarray_size_t ind1, amsarray_size_t ind2) const
+    {
+        int res;
+        amsarray<T> ret;
+        T defval = T();
+        amsarray_size_t I,J;
+
+        int nl = ind2-ind1;
+
+        if(nl<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        ret.resize(nl);
+        for(I=0;I<nl;I++)
+        {
+            J = I + ind1;
+            if(J<0 || J>=this->length)
+            {
+                ret.data[I] = defval;
+            }
+            else
+            {
+                ret.data[I] = this->data[J];
+            }
+        }
+        
+        return ret;
+    }
+
+//inserts a value at (before) index ind
+//insert(0,val) would give {val, oldval0, oldval1, ...}
+//insert(N,val)
+template<typename T> int amsarray<T>::insert(amsarray_size_t ind, const T& val)
+{
+    int ret = amsarray_success;
+    int res;
+    amsarray<T> narr;
+    amsarray_size_t I;
+
+    if(ind<0)
+    {
+        ret = amsarray_failure;
+        return ret;
+    }
+    else if(ind<=this->length)
+    {
+        res = narr.resize(this->length+1);
+        if(res!=amsarray_success)
+        {
+            ret = amsarray_failure;
+            return ret;
+        }
+
+        if(this->data!=NULL)
+        {
+            for(I=0;I<ind && I<this->length;I++)
+            {
+                narr.data[I] = this->data[I];
+            }
+        }
+        narr.data[ind] = val;
+        if(this->data!=NULL)
+        {
+            for(I=ind+1;I<narr.length;I++)
+            {
+                narr.data[I] = this->data[I-1];
+            }
+        }
+
+        //steal data from narr
+        if(this->data!=NULL) {delete[] this->data; this->data=NULL;}
+        this->length = narr.length;
+        this->data = narr.data;
+        narr.length = 0;
+        narr.data = NULL;
+    }
+    else
+    {
+        //inserting past the end of the array
+        res = narr.resize(ind+1);
+        if(res!=amsarray_success)
+        {
+            ret = amsarray_failure;
+            return ret;
+        }
+
+        if(this->data!=NULL)
+        {
+            for(I=0;I<this->length;I++)
+            {
+                narr.data[I] = this->data[I];
+            }
+        }
+        narr.data[ind] = val;
+
+        //steal data from narr
+        if(this->data!=NULL) {delete[] this->data; this->data=NULL;}
+        this->length = narr.length;
+        this->data = narr.data;
+        narr.length = 0;
+        narr.data = NULL;
+    }
+
+    return ret;
+}
+
+//removes the value at index ind, and reduces array size by 1
+template<typename T> int amsarray<T>::erase(amsarray_size_t ind)
+{
+    int ret = amsarray_success;
+    int res;
+    amsarray<T> narr;
+    amsarray_size_t I;
+
+    if(ind<0 || ind>=this->length)
+    {
+        ret = amsarray_failure;
+        return ret;
+    }
+
+    res = narr.resize(this->length-1);
+    if(res!=amsarray_success)
+    {
+        ret = amsarray_failure;
+        return ret;
+    }
+
+    if(this->data!=NULL)
+    {
+        for(I=0;I<ind && I<this->length && I<narr.length;I++)
+        {
+            narr.data[I] = this->data[I];
+        }
+        for(I=ind+1;I<this->length && (I-1)<narr.length;I++)
+        {
+            narr.data[I-1] = this->data[I];
+        }
+    }
+
+    //steal data from narr (should be similar to what std::move would do?)
+    if(this->data!=NULL) {delete[] this->data; this->data=NULL;}
+    this->length = narr.length;
+    this->data = narr.data;
+    narr.length = 0;
+    narr.data = NULL;
+
+    return ret;
+}
+
+//finds the first instance of val in the array
+template<typename T> amsarray_size_t amsarray<T>::find(const T& val)
+{
+    amsarray_size_t ret = -1;
+    amsarray_size_t I;
+
+    for(I=0;I<this->length;I++)
+    {
+        if(this->data[I]==val)
+        {
+            ret = I;
+            break;
+        }
+    }
+
+    return ret;
+}
+
+//finds the next instance of val, starting consideration at indstart
+template<typename T> amsarray_size_t amsarray<T>::findnext(amsarray_size_t indstart, const T& val)
+{
+    amsarray_size_t ret = -1;
+    amsarray_size_t I;
+
+    for(I=indstart;I<this->length;I++)
+    {
+        if(this->data[I]==val)
+        {
+            ret = I;
+            break;
+        }
+    }
+
+    return ret;
+}
+
+//finds where to insert a particular value such that the list remains ordered
+template<typename T> amsarray_size_t amsarray<T>::find_insert_ordered(const T& val)
+{
+    amsarray_size_t ret = -1;
+    amsarray_size_t I;
+
+    ret = 0;
+    for(I=0;I<this->length;I++)
+    {
+        ret = I+1;
+        if(this->data[I]>val)
+        {
+            ret = I;
+            break;
+        }
+    }
+    ret = (ret<0) ? 0 : ret;
+    ret = (ret>this->length) ? this->length : ret;
+    
+
+    return ret;
+}
+
+
+template<typename T> void amsarray_operator_comp_tf( 
+    const amsarray<T> *a,
+    const amsarray<T> *b, 
+    amsarray<bool> *cmp,
+    int threadnum, int nthreads
+)
+{
+    int I0,I1,Is,I;
+
+    
+    Is = (nthreads<=1) ? a->length : a->length/nthreads;
+    I0 = Is*threadnum;
+    I1 = (threadnum>=(nthreads-1))? a->length: Is*(threadnum+1);
+    
+    cmp->data[threadnum] = 1;
+    for(I=I0;I<I1;I++)
+    {
+        if(a->data[I]!=b->data[I])
+        {
+            cmp->data[threadnum] = 0;
+            break;
+        }
+    }
+    return;
+}
+
+template<typename T> bool amsarray<T>::operator==(const amsarray<T>& other) const
+{
+    bool ret = 1;
+    amsarray<bool> cmp;
+    amsarray_size_t I;
+    int J;
+    int nthreads;
+    std::vector<std::thread*> threads;
+
+    if(this->length != other.length)
+    {
+        ret = 0;
+        return ret;
+    }
+    if(this->length == 0)
+    {
+        ret = 1;
+        return ret;
+    }
+    if(this->length < amsmathutil25_threadpsz)
+    {
+        for(I=0;I<this->length;I++)
+        {
+            if(this->data[I]!=other.data[I])
+            {
+                ret = 0;
+                break;
+            }
+        }
+    }
+    else
+    {
+        //threaded operation
+        nthreads = std::thread::hardware_concurrency();
+        nthreads = (nthreads<1) ? 1 : nthreads;
+        nthreads = (nthreads>amsmathutil25_maxthreads) ? amsmathutil25_maxthreads : nthreads;
+        //if(nthreads<1) nthreads=1;
+        //if(nthreads>narray_max_threads) nthreads = narray_max_threads;
+        threads.resize(nthreads);
+        cmp.resize(nthreads);
+        for(J=0;J<nthreads;J++)
+        {
+            //printf("debug: launching thread %d\n",J);
+            threads[J] = new(std::nothrow) std::thread(
+                &amsarray_operator_comp_tf<T>,
+                this,
+                &other,
+                &cmp,
+                J,nthreads
+            );
+        }
+        for(J=0;J<nthreads;J++)
+        {
+            if(threads[J]!=NULL)
+            { 
+                threads[J]->join(); 
+                delete threads[J]; 
+                threads[J]= NULL;
+            }
+        }
+
+        ret = 1;
+        for(J=0;J<nthreads;J++)
+        {
+            if(cmp[J]==0)
+            {
+                ret = 0;
+                break;
+            }
+        }
+
+        cmp.resize(0);
+    }
+
+    return ret;
+}
+
+template<typename T> bool amsarray<T>::operator!=(const amsarray<T>& other) const
+{
+    return !(*this==other);
+}
+
+//aliases to insert/erase operations
+template<typename T> int amsarray<T>::append(const T& val)
+{
+    int ret = amsarray_success;
+    ret = insert(val,length);
+    return ret;
+}
+
+template<typename T> int amsarray<T>::prepend(const T& val)
+{
+    int ret = amsarray_success;
+    ret = insert(val,0);
+    return ret;
+}
+
+template<typename T> int amsarray<T>::push_back(const T& val)
+{
+    return append(val);
+}
+
+template<typename T> int amsarray<T>::push_front(const T& val)
+{
+    return prepend(val);
+}
+
+template<typename T> T amsarray<T>::pop_back()
+{
+    T ret = data[length-1];
+    erase(length-1);
+    return ret;
+}
+
+template<typename T> T amsarray<T>::pop_front()
+{
+    T ret = data[0];
+    erase(0);
+    return ret;
+}
+
+template<typename T> void amsarray_select_tf(
+    int threadnum,
+    int nthreads,
+    amsarray<T> *array,
+    amsarray<amsarray_size_t> *inds,
+    amsarray<T> *ret
+)
+{
+    amsarray_size_t I,I0,I1,Is,N;
+    amsarray_size_t ind, N2;
+    T defval = T();
+    N = inds->length;
+    N2 = array->length;
+    N = (N<=0) ? 0 : N;
+    Is = (nthreads>=1) ? N/nthreads : N;
+    I0 = Is*(threadnum);
+    I1 = (threadnum<nthreads-1) ? Is*(threadnum+1) : N;
+
+    for(I=I0;I<I1;I++)
+    {
+        ind = inds->at(I);
+        if(ind<0 || ind >=N2)
+        {
+            ret->data[I] = defval;
+        }
+        else
+        {
+            ret->data[I] = array->data[ind];
+        }
+    }
+    return;
+}
+
+//returns an array that is {thisarray[inds[0]],thisarray[inds[1]],...}
+template<typename T> amsarray<T> amsarray<T>::select(amsarray<amsarray_size_t> inds)
+{
+    int res;
+    amsarray<T> ret;
+    amsarray_size_t psize = inds.length;
+    amsarray_size_t I, ind;
+    T defval = T();
+
+    res = ret.resize(psize);
+    if(res!=amsarray_success)
+    {
+        ret.resize(0);
+        return ret;
+    }
+
+    if(psize<amsmathutil25_threadpsz)
+    {
+        for(I=0;I<psize;I++)
+        {
+            ind = inds[I];
+            if(ind<0 || ind >=length)
+            {
+                ret[I] = defval;
+            }
+            else
+            {
+                ret[I] = data[ind];
+            }
+        }
+    }
+    else
+    {
+        threaded_execute(
+            &amsarray_select_tf,
+            inds.length,
+            this,
+            &inds,
+            &ret
+        );
+    }
+
+    return ret;
+}
+
+template<typename T> void amsarray_reverse_tf(
+    int threadnum,
+    int nthreads,
+    amsarray<T> *array,
+    amsarray<T> *ret
+)
+{
+    amsarray_size_t I,I0,I1,Is,N;
+    T defval = T();
+    N = array->length;
+    N = (N<=0) ? 0 : N;
+    Is = (nthreads>=1) ? N/nthreads : N;
+    I0 = Is*(threadnum);
+    I1 = (threadnum<nthreads-1) ? Is*(threadnum+1) : N;
+
+    for(I=I0;I<I1;I++)
+    {
+        ret->data[I] = array->data[N-I-1];
+    }
+
+    return;
+}
+
+//returns an array that is this array in reverse order
+template<typename T> amsarray<T> amsarray<T>::reverse()
+{
+    int res;
+    amsarray<T> ret;
+    amsarray_size_t psize = length;
+    amsarray_size_t I;
+    T defval = T();
+
+    res = ret.resize(psize);
+    if(res!=amsarray_success)
+    {
+        ret.resize(0);
+        return ret;
+    }
+
+    if(psize<amsmathutil25_threadpsz)
+    {
+        for(I=0;I<psize;I++)
+        {
+            ret[I] = data[length-I-1];
+        }
+    }
+    else
+    {
+        threaded_execute(
+            &amsarray_reverse_tf,
+            length,
+            this,
+            &ret
+        );
+    }
+
+    return ret;
+}
+
+
+
+
+
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/util/amsmathutil25_amsarray_sortimpl.hpp b/include/amsmathutil25/util/amsmathutil25_amsarray_sortimpl.hpp
new file mode 100644
index 0000000..a5e68a9
--- /dev/null
+++ b/include/amsmathutil25/util/amsmathutil25_amsarray_sortimpl.hpp
@@ -0,0 +1,90 @@
+#ifndef __AMSMATHUTIL25_AMSARRAY_SORTIMPL_HPP__
+#define __AMSMATHUTIL25_AMSARRAY_SORTIMPL_HPP__
+
+namespace ams
+{
+
+
+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,
+    ams::pair<amsarray_size_t,amsarray_size_t> *leftrange,
+    ams::pair<amsarray_size_t,amsarray_size_t> *rightrange
+)
+{
+    int ret = 0;
+    bool b1,b2;
+    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)
+    {
+        //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;
+        return ret;
+    }
+    if((range.b - range.a) == 2)
+    {
+        //two element range - sort directly
+        v1 = array->data[permarray->data[range.a]];
+        v2 = array->data[permarray->data[range.b]];
+        if(v2<v1)
+        {
+            //swap permutation indices
+            tmp = permarray->data[range.a];
+            permarray->data[range.a] = permarray->data[range.b];
+            permarray->data[range.b] = 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;
+        return ret;
+    }
+    else
+    {
+        //perform quicksort iteration
+
+        //choose midpoint pivot
+        P = (range.a + range.b)/2;
+        P = (P<range.a) ? range.a : P;
+        P = (P>=range.b) ? range.b-1 : P;
+
+        
+
+    }
+
+
+
+
+    return ret;
+}
+
+
+//returns an array of indices that is a permutation which will sort
+//this array in ascending order
+template<typename T> amsarray<amsarray_size_t> amsarray<T>::sort_permutation()
+{
+    int res;
+    amsarray<amsarray_size_t> ret;
+
+    return ret;
+}
+
+
+
+
+
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/util/amsmathutil25_bufferops.hpp b/include/amsmathutil25/util/amsmathutil25_bufferops.hpp
new file mode 100644
index 0000000..d272583
--- /dev/null
+++ b/include/amsmathutil25/util/amsmathutil25_bufferops.hpp
@@ -0,0 +1,26 @@
+#ifndef __AMSMATHUTIL25_BUFFEROPS_HPP__
+#define __AMSMATHUTIL25_BUFFEROPS_HPP__
+
+namespace ams
+{
+    typedef int64_t buffer_size_t;
+
+    //threaded set - sets N elements of bufffer to val
+    template<typename T> int buffer_set(T* buffer, buffer_size_t N, const T val);
+
+    //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 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);
+
+    //threaded copy of bufferfrom to bufferto
+    template<typename T1, typename T2> int buffer_cast_copy(T1* bufferto, const T2* bufferfrom, buffer_size_t N);
+
+}; //end namespace ams
+
+#include <amsmathutil25/util/amsmathutil25_bufferops_impl.hpp>
+
+#endif
+
diff --git a/include/amsmathutil25/util/amsmathutil25_bufferops_impl.hpp b/include/amsmathutil25/util/amsmathutil25_bufferops_impl.hpp
new file mode 100644
index 0000000..4768370
--- /dev/null
+++ b/include/amsmathutil25/util/amsmathutil25_bufferops_impl.hpp
@@ -0,0 +1,161 @@
+#ifndef __AMSMATHUTIL25_BUFFEROPS_IMPL_HPP__
+#define __AMSMATHUTIL25_BUFFEROPS_IMPL_HPP__
+
+namespace ams
+{
+
+    template<typename T> void buffer_set_tf(
+        int threadnum,
+        int nthread,
+        T *buffer,
+        buffer_size_t indstart, 
+        buffer_size_t indstop, 
+        const T val
+    )
+    {
+        buffer_size_t I,I0,I1,Is,N;
+
+        N = indstop-indstart;
+        N = (N<0) ? 0 : N;
+        //Is = N/nthread;
+
+        Is = (nthread<=0) ? N : N/nthread;
+        I0 = Is*(threadnum);
+        I1 = (threadnum>=(nthread-1)) ? N : Is*(threadnum+1);
+        
+        // I0 = (I0<=0) ? 0 : I0;
+        // I1 = (I1<=0) ? 0 : I1;
+        I0 = (I0>N) ? N : I0;
+        I1 = (I1>N) ? N : I1;
+
+        for(I=I0;I<I1;I++)
+        {
+            buffer[I+indstart] = val;
+        }
+        
+        return;
+    }
+
+    //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)
+    {
+        int ret = amsmathutil25_success;
+        buffer_size_t psize;
+        buffer_size_t I;
+
+        indstart = (indstart<0) ? 0 : indstart;
+        indstop = (indstop<0) ? 0 : indstop;
+        indstop = (indstop<indstart) ? indstart : indstop;
+        psize = indstop-indstart;
+        //psize = (psize<=0) ? 0 : psize; //redundant
+
+        if(psize<amsmathutil25_threadpsz)
+        {
+            for(I=0;I<psize;I++)
+            {
+                buffer[I+indstart] = val;
+            }
+        }
+        else
+        {
+
+            ams::threaded_execute(
+                &buffer_set_tf<T>,(int64_t) psize,
+                buffer,indstart,indstop,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)
+    {
+        int ret = amsmathutil25_success;
+        buffer_size_t I;
+
+        if(N<amsmathutil25_threadpsz)
+        {
+            for(I=0;I<N;I++) buffer[I] = val;
+        }
+        else
+        {
+            ret = buffer_set<T>(buffer,0,N,val);
+        }
+        return ret;
+    }
+
+    template<typename T1, typename T2> void buffer_cast_copy_tf(
+        int threadnum,
+        int nthread,
+        T1 *bufferto,
+        const T2 *bufferfrom,
+        buffer_size_t offsetto,
+        buffer_size_t offsetfrom,
+        buffer_size_t N
+    )
+    {
+        buffer_size_t I,I0,I1,Is;
+
+        // Is = N/(nthread-1);
+        // Is = (Is<=0) ? 1 : Is;
+        // I0 = Is*(threadnum);
+        // I1 = Is*(threadnum+1);
+        // I0 = (I0<=0) ? 0 : I0;
+        // I1 = (I1<=0) ? 0 : I1;
+        // I0 = (I0>N) ? N : I0;
+        // I1 = (I1>N) ? N : I1;
+
+        Is = (nthread<=0) ? N : N/nthread;
+        I0 = Is*(threadnum);
+        I1 = (threadnum>=(nthread-1)) ? N : Is*(threadnum+1);
+        
+        // I0 = (I0<=0) ? 0 : I0;
+        // I1 = (I1<=0) ? 0 : I1;
+        I0 = (I0>N) ? N : I0;
+        I1 = (I1>N) ? N : I1;
+
+        for(I=I0;I<I1;I++)
+        {
+            bufferto[I + offsetto] = (T1) (bufferfrom[I + offsetfrom]);
+        }
+        
+        return;
+    }
+
+    //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)
+    {
+        int ret = amsmathutil25_success;
+
+        ams::threaded_execute(
+            &buffer_cast_copy_tf<T1,T2>,(int64_t) N,
+            bufferto,bufferfrom,offsetto,offsetfrom,N
+        );
+
+        return ret;
+    }
+
+    //threaded copy of bufferfrom to bufferto
+    template<typename T1, typename T2> int buffer_cast_copy(T1* bufferto, const T2* bufferfrom, buffer_size_t N)
+    {
+        int ret = amsmathutil25_success;
+        buffer_size_t I;
+
+        if(N<amsmathutil25_threadpsz)
+        {
+            for(I=0;I<N;I++) bufferto[I] = (T1) bufferfrom[I];
+        }
+        else
+        {
+            ret = buffer_cast_copy<T1,T2>(bufferto,bufferfrom,0,0,N);
+        }
+
+        return ret;
+    }
+
+}; //end namespace ams
+
+#endif
+
diff --git a/include/amsmathutil25/util/amsmathutil25_util.hpp b/include/amsmathutil25/util/amsmathutil25_util.hpp
new file mode 100644
index 0000000..3041858
--- /dev/null
+++ b/include/amsmathutil25/util/amsmathutil25_util.hpp
@@ -0,0 +1,38 @@
+#ifndef __AMSMATHUTIL25_UTIL_HPP__
+#define __AMSMATHUTIL25_UTIL_HPP__
+
+namespace ams
+{
+
+//A template function that takes as input a function pointer and a series of arguments
+//The function is executed with fptr(threadnum, nthreads, otherargs...) with a certain number of threads
+//psize must be supplied to determine whether to execute in threaded mode or not.
+template<typename callable, typename ... argst> int threaded_execute(callable &&fptr, int64_t psize, argst&&... args);
+
+template<typename T1, typename T2> struct pair
+{
+public:
+    T1 a;
+    T2 b;
+    pair() {a = T1(); b = T2();}
+    pair(const T1 &_a, const T2& _b) {a = _a; b = _b;}
+};
+
+template<typename T1, typename T2, typename T3> struct triple
+{
+public:
+    T1 a;
+    T2 b;
+    T3 c;
+    triple() {a = T1(); b = T2(); c = T3();}
+    triple(const T1 &_a, const T2& _b, const T3& _c) {a = _a; b = _b; c = _c;}
+};
+
+}; //end namespace ams
+
+#include <amsmathutil25/util/amsmathutil25_utilimpl.hpp>
+#include <amsmathutil25/util/amsmathutil25_bufferops.hpp>
+#include <amsmathutil25/util/amsmathutil25_amsarray.hpp>
+
+#endif
+
diff --git a/include/amsmathutil25/util/amsmathutil25_utilimpl.hpp b/include/amsmathutil25/util/amsmathutil25_utilimpl.hpp
new file mode 100644
index 0000000..0bce248
--- /dev/null
+++ b/include/amsmathutil25/util/amsmathutil25_utilimpl.hpp
@@ -0,0 +1,71 @@
+#ifndef __AMSMATHUTIL25_UTILIMPL_HPP__
+#define __AMSMATHUTIL25_UTILIMPL_HPP__
+
+namespace ams
+{
+
+//A template function that takes as input a function pointer and a series of arguments
+//The function is executed with fptr(threadnum, nthreads, otherargs...) with a certain number of threads
+//psize must be supplied to determine whether to execute in threaded mode or not.
+template<typename callable, typename ... argst> int threaded_execute(callable &&fptr, int64_t psize, argst&&... args)
+{
+    int ret = amsmathutil25_success;
+    int nthreads;
+    int I;
+    std::vector<std::thread*> threads;
+
+    if(psize<amsmathutil25_threadpsz)
+    {
+        nthreads = 1;
+        I = 0;
+        std::invoke(
+            std::forward<callable>(fptr),
+            I,
+            nthreads,
+            std::forward<argst>(args)...
+        );
+    }
+    else
+    {
+        nthreads = std::thread::hardware_concurrency(); //number of cpu cores the system thinks you have
+        nthreads = (nthreads<=0) ? 1 : nthreads;
+        nthreads = (nthreads>amsmathutil25_maxthreads) ? amsmathutil25_maxthreads : nthreads;
+        threads.resize(nthreads);
+        for(I=0;I<nthreads;I++) threads[I] = NULL;
+        for(I=0;I<nthreads;I++)
+        {
+            threads[I] = new(std::nothrow) std::thread
+            (
+                std::forward<callable>(fptr),
+                I,
+                nthreads,
+                std::forward<argst>(args)...
+            );
+        }
+        for(I=0;I<nthreads;I++)
+        {
+            if(threads[I]==NULL)
+            { //null thread creation failure check
+                //printf("debug check!\n");
+                ret = amsmathutil25_failure;
+            }
+        }
+        for(I=0;I<nthreads;I++)
+        {
+            if(threads[I]!=NULL)
+            {
+                threads[I]->join();
+                delete threads[I];
+                threads[I] = NULL;
+            }
+        }
+    }
+
+    return ret;
+}
+
+}; //end namespace ams
+
+
+#endif
+
diff --git a/src/amsmathutil25/testing/amsmathtuil25_test1.cpp b/src/amsmathutil25/testing/amsmathtuil25_test1.cpp
new file mode 100644
index 0000000..696873f
--- /dev/null
+++ b/src/amsmathutil25/testing/amsmathtuil25_test1.cpp
@@ -0,0 +1,96 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+namespace amsmathutil25
+{
+
+void test_amsarray1()
+{
+    amsarray<float> a,b,c,q;
+    int64_t I;
+
+    printf("Tests for amsarray...\n");
+    q.resize(10);
+    q[0] = 1;
+    q[9] = 5;
+    for(I=0;I<10;I++)
+    {
+        printf("q[%ld] = %1.3f\n",I,q[I]);
+    }
+
+    q.resize(12);
+    a = q;
+    a = b = q;
+    for(I=0;I<12;I++)
+    {
+        printf("a[%ld] = %1.3f\n",I,a[I]);
+    }
+
+    q.resize(8);
+    a = q;
+    a = b = q;
+    for(I=0;I<8;I++)
+    {
+        printf("a[%ld] = %1.3f\n",I,a[I]);
+    }
+    
+    a.resize(amsmathutil25_threadpsz + 10);
+    a.setall(55);
+    printf("a.length=%ld\n",a.length);
+
+    b = a;
+    for(I=amsmathutil25_threadpsz;I<b.length;I++)
+    {
+        printf("b[%ld] = %1.3f\n",I,b[I]);
+    }
+
+    return;
+}
+
+void test_amsarray2()
+{
+    amsarray<int64_t> a,b,c;
+    amsarray<float> q;
+    int64_t I;
+
+    printf("Tests for amsarray...\n");
+
+    a.resize(0);
+    b = a;
+    a.resize(-10);
+    b = a;
+
+    printf("a==b?%d\n",a==b);
+
+    a.resize(10);
+    b = a;
+
+    printf("a==b?%d\n",a==b);
+
+    a.resize(10000);
+    b = a;
+
+    printf("a==b?%d\n",a==b);
+
+    a.setall(55);
+    b = a;
+
+    printf("a==b?%d\n",a==b);
+
+    a[0] = 1;
+
+    printf("a==b?%d\n",a==b);
+
+    q = (amsarray<float>)a;
+    for(I=q.length-10;I<q.length;I++)
+    {
+        printf("q[%ld]=%1.3f\n",I,q[I]);
+    }
+
+
+}
+
+
+}; //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
new file mode 100644
index 0000000..663a136
--- /dev/null
+++ b/src/amsmathutil25/util/amsmathutil25_amsarray.cpp
@@ -0,0 +1,84 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+
+template<> void amsarray<int>::print(bool newline,int printstyle)
+{
+    amsarray_size_t I;
+
+    printf("{");
+    if(data!=NULL)
+    {
+        for(I=0;I<length-1;I++)
+        {
+            printf("%d,",data[I]);
+        }
+        if(length>0) printf("%d",data[length-1]);
+    }
+    printf("}");
+    if(newline==1) printf("\n");
+
+    return;
+}
+
+template<> void amsarray<long>::print(bool newline,int printstyle)
+{
+    amsarray_size_t I;
+
+    printf("{");
+    if(data!=NULL)
+    {
+        for(I=0;I<length-1;I++)
+        {
+            printf("%ld,",data[I]);
+        }
+        if(length>0) printf("%ld",data[length-1]);
+    }
+    printf("}");
+    if(newline==1) printf("\n");
+
+    return;
+}
+
+template<> void amsarray<float>::print(bool newline,int printstyle)
+{
+    amsarray_size_t I;
+
+    printf("{");
+    if(data!=NULL)
+    {
+        for(I=0;I<length-1;I++)
+        {
+            printf("%1.3f,",data[I]);
+        }
+        if(length>0) printf("%1.3f",data[length-1]);
+    }
+    printf("}");
+    if(newline==1) printf("\n");
+
+    return;
+}
+
+template<> void amsarray<double>::print(bool newline,int printstyle)
+{
+    amsarray_size_t I;
+
+    printf("{");
+    if(data!=NULL)
+    {
+        for(I=0;I<length-1;I++)
+        {
+            printf("%1.3f,",data[I]);
+        }
+        if(length>0) printf("%1.3f",data[length-1]);
+    }
+    printf("}");
+    if(newline==1) printf("\n");
+
+    return;
+}
+
+
+
+};
\ No newline at end of file
diff --git a/src/amsmathutil25/util/amsmathutil25_amsarray_sort.cpp b/src/amsmathutil25/util/amsmathutil25_amsarray_sort.cpp
new file mode 100644
index 0000000..8b7af57
--- /dev/null
+++ b/src/amsmathutil25/util/amsmathutil25_amsarray_sort.cpp
@@ -0,0 +1,62 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+namespace ams
+{
+
+void amsarray_permutation_identity_tf(
+    int threadnum,
+    int nthreads,
+    amsarray<amsarray_size_t> *ret
+)
+{
+    amsarray_size_t I,I0,I1,Is,N;
+    N = ret->length;
+    N = (N<=0) ? 0 : N;
+    Is = (nthreads>=1) ? N/nthreads : N;
+    I0 = Is*(threadnum);
+    I1 = (threadnum<nthreads-1) ? Is*(threadnum+1) : N;
+
+    for(I=I0;I<I1;I++)
+    {
+        ret->data[I] = I;
+    }
+
+
+    return;
+}
+
+// {0, 1, 2, .... N}
+amsarray<amsarray_size_t> permutation_identity(amsarray_size_t _length)
+{
+    amsarray<amsarray_size_t> ret;
+    int res;
+    amsarray_size_t psize = _length;
+    amsarray_size_t I;
+    
+    res = ret.resize(psize);
+    if(res!=amsarray_success)
+    {
+        ret.resize(0);
+        return ret;
+    }
+
+    if(psize<amsmathutil25_threadpsz)
+    {
+        for(I=0;I<psize;I++)
+        {
+            ret[I] = I;
+        }
+    }
+    else
+    {
+        threaded_execute(
+            &amsarray_permutation_identity_tf,
+            _length,
+            &ret
+        );
+    }
+
+    return ret;
+}
+
+};
\ No newline at end of file
diff --git a/src/main.cpp b/src/main.cpp
index 46654d6..97ae37e 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -5,4 +5,7 @@ int main(int argc, char* argv[])
     int ret = 0;
     printf("ams c++ math and utility library tests.\n");
 
+    //ams::amsmathutil25::test_amsarray1();
+    //ams::amsmathutil25::test_amsarray2();
+    return ret;
 }
\ No newline at end of file