random updates

amsmathutil25.code-workspace

@@ -8,6 +8,9 @@
 		},
 		{
 			"path": "../amscppnarray"
+		},
+		{
+			"path": "../../sourceprojs23/amscpptimelib"
 		}
 	]
 }

include/amsmathutil25/math/amsmathutil25_complex128.hpp

@@ -13,7 +13,7 @@ namespace cmp
 
         complex();
         complex(double _r, double _i);
-        complex(double _r);
+        explicit complex(double _r);
         complex(const complex &other);
         complex& operator=(const complex& other);
         complex& operator=(const double other);

include/amsmathutil25/math/amsmathutil25_complex64.hpp

@@ -13,7 +13,7 @@ namespace cmp
 
         complex64();
         complex64(float _r, float _i);
-        complex64(float _r);
+        explicit complex64(float _r);
         complex64(const complex64 &other);
         complex64& operator=(const complex64& other);
         complex64& operator=(const float other);

include/amsmathutil25/random/amsmathutil25_random.hpp

@@ -11,17 +11,36 @@ 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_randt1 dpr32_add1 = ( (amsmu_randt1) 293487 );
 
 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 );
+static const amsmu_randt2 dpr64_add1 = ( (amsmu_randt2) 293487L );
 
 extern amsmu_randt1 dpr32_rseed;
 extern amsmu_randt2 dpr64_rseed;
 
-amsmu_randt1 dpr32_nextseed(amsmu_randt1 seed);
-amsmu_randt2 dpr64_nextseed(amsmu_randt2 seed);
+amsmu_randt1 seed32_next(amsmu_randt1 seed);
+amsmu_randt2 seed64_next(amsmu_randt2 seed);
+
+double rand(amsmu_randt1 *seed = &dpr32_rseed);
+float randf(amsmu_randt1 *seed = &dpr32_rseed);
+double randgaussian(amsmu_randt1 *seed = &dpr32_rseed);
+float randgaussianf(amsmu_randt1 *seed = &dpr32_rseed);
+int randint(int low, int highexcl, amsmu_randt1 *seed = &dpr32_rseed);
+int64_t randintl(int64_t low, int64_t highexcl, amsmu_randt2 *seed = &dpr64_rseed);
+
+void seed32_set(amsmu_randt1 _seed);
+void seed64_set(amsmu_randt2 _seed);
+void seed_init_timer();
+
+//Threaded generation of random amsarrays
+
+amsarray<double> amsarray_rand(amsarray_size_t N, amsmu_randt1 *seed = &dpr32_rseed);
+amsarray<float> amsarray_randf(amsarray_size_t N, amsmu_randt1 *seed = &dpr32_rseed);
+
+amsarray<int> amsarray_randint(amsarray_size_t N, int low, int highexcl, amsmu_randt1 *seed = &dpr32_rseed);
+amsarray<int64_t> amsarray_randintl(amsarray_size_t N, int64_t low, int64_t highexcl, amsmu_randt2 *seed = &dpr64_rseed);
 
 
 

include/amsmathutil25/testing/amsmathutil25_testing.hpp

@@ -8,9 +8,11 @@ namespace amsmathutil25
 
     void test_amsarray1();
     void test_amsarray2();
-
-    void test_select();
     void test_amsarray_sort();
+    void test_amsarray_select();
+    void test_amsarray_sort1();
+
+    
 
 }; //end namespace amsmathutil25
 }; //end namespace ams

include/amsmathutil25/util/amsmathutil25_amsarray.hpp

@@ -7,6 +7,8 @@ namespace ams
     static const int amsarray_success = amsmathutil25_success;
     static const int amsarray_failure = amsmathutil25_failure;
 
+    static const int amsarray_sortthreadpsize = 5000;
+
     template<typename T> class amsarray
     {
         public:
@@ -109,6 +111,9 @@ namespace ams
 // {0, 1, 2, .... N}
 amsarray<amsarray_size_t> permutation_identity(amsarray_size_t _length);
 
+void _debug_amsarray_print(amsarray<double> *array, bool newline=1,int printstyle=0);
+void _debug_amsarray_print(amsarray<long> *array, bool newline=1,int printstyle=0);
+
 
 }; //end namespace ams
 

include/amsmathutil25/util/amsmathutil25_amsarray_sortimpl.hpp

@@ -41,9 +41,10 @@ template<typename T> int amsarray_quicksort_round(
         if(v2<v1)
         {
             //swap permutation indices
-            tmp = permarray->data[range.a];
-            permarray->data[range.a] = permarray->data[range.b-1];
-            permarray->data[range.b-1] = tmp;
+            amsarray_permutation_swap(permarray,range.a,range.b-1);
+            // tmp = permarray->data[range.a];
+            // 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
         *leftrange = ams::pair<amsarray_size_t,amsarray_size_t>(-1,-1);
@@ -63,18 +64,28 @@ template<typename T> int amsarray_quicksort_round(
         amsarray_permutation_swap(permarray,P,range.b-1);
         P = range.b-1;
         
+        // amsarray_permutation_swap(permarray,P,range.a);
+        // P = range.a;
+        
+
         J = range.a;
         for(I=range.a;I<range.b-1;I++)
         {
+            
+            //printf("debug: I=%ld, J=%ld, P=%ld, a=%1.3f, b=%1.3f\n",I,J,P,
+            //    (double)array->data[permarray->data[I]],(double)array->data[permarray->data[P]]);
+
             if(array->data[permarray->data[I]]<array->data[permarray->data[P]])
             {
                 if(J!=I)
                 {
+                    //printf("debug: swap\n");
                     amsarray_permutation_swap(permarray,I,J);
                     J++;
                 }
                 else
                 {
+                    //printf("debug: skip\n");
                     J++;
                 }
             }
@@ -82,7 +93,7 @@ template<typename T> int amsarray_quicksort_round(
 
         if(array->data[permarray->data[J]]<array->data[permarray->data[P]])
         {
-            J++;
+            J++;  
             amsarray_permutation_swap(permarray,P,J);
 
         }
@@ -119,6 +130,42 @@ template<typename T> int amsarray_quicksort_round(
     return ret;
 }
 
+template<typename T> int amsarray_quicksort_subrange(
+    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
+)
+{
+    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;
+
+    ranges.append(_range);
+    rangeptr = 0;
+    while(rangeptr<ranges.length)
+    {
+        //printf("debug2:"); _debug_amsarray_print(permarray);
+        range = ranges[rangeptr];
+        rangeptr++;
+        //printf("debug3: range=(%d,%d)\n",(int)range.a,(int)range.b);
+        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> int amsarray_quicksort_unthreaded(
     amsarray<T> *array,                   //size N - array to sort
     amsarray<amsarray_size_t> *permarray //size N - permutation of sorting
@@ -145,8 +192,10 @@ template<typename T> int amsarray_quicksort_unthreaded(
     rangeptr = 0;
     while(rangeptr<ranges.length)
     {
+        //printf("debug2:"); _debug_amsarray_print(permarray);
         range = ranges[rangeptr];
         rangeptr++;
+        //printf("debug3: range=(%d,%d)\n",(int)range.a,(int)range.b);
         amsarray_quicksort_round(array,permarray,range,&rangeleft,&rangeright);
         if(rangeleft.a>=0 && rangeleft.b>rangeleft.a)
         {
@@ -173,31 +222,48 @@ template<typename T> void amsarray_quicksort_tf(
 
     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);
-        }
+    if(range.b-range.a < amsarray_sortthreadpsize)
+    {
+        res = amsarray_quicksort_subrange(
+            array,
+            permarray,
+            range
+        );
+        //there should be no work to be done after this
+    }
+    else
+    {
+        //range is too big, quicksort the pivot and supply subranges
+        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)
     }
-    //end critical section (end of function)
 
     return;
 }
 
+//TODO - if the range falls below a specified size, I want to be able to run through
+// quicksorting the entire range within a thread before returning
+
+
 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
@@ -240,9 +306,13 @@ template<typename T> int amsarray_quicksort_threaded(
         nthreads = ranges.length-rangeptr;
         nthreads = (nthreads>maxthreads) ? maxthreads : nthreads;
 
+        //printf("debug: %d %d %ld %ld\n",nthreads,maxthreads,rangeptr,ranges.length);
         for(I=0;I<nthreads;I++)
         {
-            range = ranges[rangeptr+I];
+            threadlock.lock();
+            range = ranges[rangeptr]; rangeptr++;
+            threadlock.unlock();
+            //printf("debug: thread %ld exec with range(%ld,%ld), rptr=%ld rlen=%ld\n",I,range.a,range.b,rangeptr,ranges.length);
             threads[I] = new(std::nothrow) std::thread(
                 amsarray_quicksort_tf<T>,
                 array,permarray,
@@ -259,6 +329,7 @@ template<typename T> int amsarray_quicksort_threaded(
                 ret = amsarray_failure;
             }
         }
+        //printf("debug3\n");
         for(I=0;I<nthreads;I++)
         {
             if(threads[I]!=NULL)
@@ -267,7 +338,7 @@ template<typename T> int amsarray_quicksort_threaded(
                 delete threads[I];
                 threads[I] = NULL;
             }
-            rangeptr++;
+            //rangeptr++;
         }
     }
 
@@ -292,7 +363,7 @@ template<typename T> int amsarray_quicksort(
         }
     }
 
-    if(array->length<amsmathutil25_threadpsz)
+    if(array->length<amsarray_sortthreadpsize)
     {
         //perform unthreaded quicksort
         ret = amsarray_quicksort_unthreaded(

include/amsmathutil25/util/amsmathutil25_util.hpp

@@ -28,6 +28,9 @@ public:
     triple(const T1 &_a, const T2& _b, const T3& _c) {a = _a; b = _b; c = _c;}
 };
 
+//returns time in msec
+double time_msec();
+
 }; //end namespace ams
 
 #include <amsmathutil25/util/amsmathutil25_utilimpl.hpp>

src/amsmathutil25/math/amsmathutil25_mathfns1.cpp

@@ -2,8 +2,6 @@
 
 namespace ams
 {
-namespace amsmathutil25
-{
 
 bool isnan(double d)
 {
@@ -210,5 +208,4 @@ int64_t abs(int64_t a)
     return ret;
 }
 
-};
-};
+}; //end namespace ams

src/amsmathutil25/random/amsmathutiil25_random.cpp

@@ -8,22 +8,387 @@ namespace rand
 amsmu_randt1 dpr32_rseed;
 amsmu_randt2 dpr64_rseed;
 
-amsmu_randt1 dpr32_nextseed(amsmu_randt1 seed)
+amsmu_randt1 seed32_next(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 seed64_next(amsmu_randt2 seed)
 {
     amsmu_randt2 sret = seed;
     sret = ams::mod(sret*dpr64_mult1+dpr64_add1,dpr64_mod);
     return sret;
 }
 
+double rand(amsmu_randt1 *seed)
+{
+    double ret = 0.0;
+    *seed = seed32_next(*seed);
+    ret = ((double) *seed) / ((double)(dpr32_mod-1));
+    return ret;
+}
+
+float randf(amsmu_randt1 *seed)
+{
+    float ret = 0.0;
+    *seed = seed32_next(*seed);
+    ret = ((float) *seed) / ((float)(dpr32_mod-1));
+    return ret;
+}
+
+double randgaussian(amsmu_randt1 *seed)
+{
+    double ret = 0.0;
+    double u1,u2;
+    u1 = rand(seed);
+    u2 = rand(seed);
+    if(u1>0.0)
+    {
+        ret = ::sqrt(-2.0*::log(u1))*::cos(2.0*ams::pi*u2);
+    }
+    return ret;
+}
+
+float randgaussianf(amsmu_randt1 *seed)
+{
+    float ret = 0.0;
+    float u1,u2;
+    u1 = randf(seed);
+    u2 = randf(seed);
+    if(u1>0.0f)
+    {
+        ret = ::sqrtf(-2.0*::logf(u1))*::cosf(2.0f*ams::pif*u2);
+    }
+    return ret;
+}
+
+int randint(int low, int highexcl, amsmu_randt1 *seed)
+{
+    int ret = 0;
+
+    if(highexcl-low>0)
+    {
+        *seed = seed32_next(*seed);
+        ret = low + (int)((*seed)%(highexcl-low));
+    }
+
+    return ret;
+}
+
+int64_t randintl(int64_t low, int64_t highexcl, amsmu_randt2 *seed)
+{
+    int64_t ret = 0;
+
+    if(highexcl-low>0)
+    {
+        *seed = seed64_next(*seed);
+        ret = low + (int)((*seed)%(highexcl-low));
+    }
+
+    return ret;
+}
+
+void seed32_set(amsmu_randt1 _seed)
+{
+    dpr32_rseed = _seed;
+    return;
+}
+
+void seed64_set(amsmu_randt2 _seed)
+{
+    dpr64_rseed = _seed;
+    return;
+}
+
+void seed_init_timer()
+{
+    amsmu_randt1 t1 = (amsmu_randt1)time(NULL);
+    amsmu_randt1 t2 = (amsmu_randt1)(::fmod((double)clock()/((double)CLOCKS_PER_SEC)*1000.0f,36000.0f));
+    dpr32_rseed = (amsmu_randt1)t1 + (amsmu_randt1)t2;
+
+    amsmu_randt2 t3 = (amsmu_randt2)time(NULL);
+    amsmu_randt2 t4 = (amsmu_randt2)(::fmod((double)clock()/((double)CLOCKS_PER_SEC)*1000.0f,36000.0f));
+    dpr64_rseed = (amsmu_randt2)t3 + (amsmu_randt2)t4;
+
+    return;
+}
+
+//Threaded generation of random amsarrays
+
+template<typename T> void amsarray_rand_threadf1(
+    amsarray<T> *out,
+    amsarray<amsmu_randt1> *rseeds,
+    T (*randfunc)(amsmu_randt1 *),
+    int threadnum,
+    int nthreads
+)
+{
+    int I0,I1,Is,I;
+    Is = out->length/nthreads;
+    I0 = Is*threadnum;
+    I1 = (threadnum>=(nthreads-1)) ? out->length : Is*(threadnum+1);
+    for(I=I0;I<I1;I++)
+    {
+        out->data[I] = randfunc(&(rseeds->data[threadnum]));
+    }
+    return;
+}
+
+template<typename T> void amsarray_random_threadexec1(
+    amsarray<T> *out,
+    amsarray_size_t N,
+    T (*randfunc)(amsmu_randt1 *),
+    amsmu_randt1 *rseed
+)
+{
+    amsarray_size_t I;
+    int J;
+    int nthreads;
+    std::vector<std::thread*> threads;
+    ams::amsarray<amsmu_randt1> rseeds;
+    int res;
+
+    res = out->resize(N);
+    if(res!=amsarray_success)
+    {
+        out->resize(0);
+        return;
+    }
+
+    if(N<amsmathutil25_threadpsz)
+    {
+        //single threaded
+        for(I=0;I<N;I++)
+        {
+            out->data[I] = randfunc(rseed);
+        }
+    }
+    else
+    {
+        //multi-threaded operation
+        nthreads = std::thread::hardware_concurrency();
+        nthreads = (nthreads<1) ? 1 : nthreads;
+        nthreads = (nthreads>amsmathutil25_maxthreads) ? amsmathutil25_maxthreads : nthreads;
+        threads.resize(nthreads);
+        rseeds.resize(nthreads);
+        
+        for(J=0;J<nthreads;J++)
+        {
+            *rseed = seed32_next(*rseed);
+            rseeds.data[J] = *rseed + 13*J;
+        }
+        for(J=0;J<nthreads;J++)
+        {
+            threads[J] = new(std::nothrow) std::thread(
+                &amsarray_rand_threadf1<T>,
+                out,
+                &rseeds,
+                randfunc,
+                J,nthreads
+            );
+            if(threads[J]==NULL)
+            {
+                //handle thread allocation error
+                printf("warning: amsarray_random_threadexec1:: thread %d failed to allocate.\n",J);
+            }
+        }
+        for(J=0;J<nthreads;J++)
+        {
+            if(threads[J]!=NULL)
+            { 
+                threads[J]->join(); 
+                delete threads[J]; 
+                threads[J]= NULL;
+            }
+        }
+    }
+
+    return;
+}
 
 
+amsarray<double> amsarray_rand(amsarray_size_t N, amsmu_randt1 *seed)
+{
+    amsarray<double> ret;
 
-};
-};
+    if(N<=0)
+    {
+        ret.resize(0);
+        return ret;
+    }
+
+    amsarray_random_threadexec1(
+        &ret,
+        N,
+        &(rand),
+        seed
+    );
+
+    return ret;
+}
+
+amsarray<float> amsarray_randf(amsarray_size_t N, amsmu_randt1 *seed)
+{
+    amsarray<float> ret;
+
+    if(N<=0)
+    {
+        ret.resize(0);
+        return ret;
+    }
+
+    amsarray_random_threadexec1(
+        &ret,
+        N,
+        &(randf),
+        seed
+    );
+
+    return ret;
+}
+
+template<typename T, typename T2> void amsarray_rand_threadf2(
+    amsarray<T> *out,
+    amsarray<T2> *rseeds,
+    T (*randfunc)(T, T, T2*),
+    T low,
+    T highexcl,
+    int threadnum,
+    int nthreads
+)
+{
+    int I0,I1,Is,I;
+    Is = out->length/nthreads;
+    I0 = Is*threadnum;
+    I1 = (threadnum>=(nthreads-1)) ? out->length : Is*(threadnum+1);
+    for(I=I0;I<I1;I++)
+    {
+        out->data[I] = randfunc(low,highexcl,&(rseeds->data[threadnum]));
+    }
+    return;
+}
+
+template<typename T, typename T2> void amsarray_random_threadexec2(
+    amsarray<T> *out,
+    amsarray_size_t N,
+    T (*randfunc)(T, T, T2*),
+    T2 (*nextfunc)(T2),
+    T low,
+    T highexcl,
+    T2 *rseed
+)
+{
+    amsarray_size_t I;
+    int J;
+    int nthreads;
+    std::vector<std::thread*> threads;
+    amsarray<T2> rseeds;
+    int res;
+
+    res = out->resize(N);
+    if(res!=amsarray_success)
+    {
+        out->resize(0);
+        return;
+    }
+
+    if(N<amsmathutil25_threadpsz)
+    {
+        //single threaded
+        for(I=0;I<N;I++)
+        {
+            out->data[I] = randfunc(low,highexcl,rseed);
+        }
+    }
+    else
+    {
+        //multi-threaded operation
+        nthreads = std::thread::hardware_concurrency();
+        nthreads = (nthreads<1) ? 1 : nthreads;
+        nthreads = (nthreads>amsmathutil25_maxthreads) ? amsmathutil25_maxthreads : nthreads;
+        threads.resize(nthreads);
+        rseeds.resize(nthreads);
+        for(J=0;J<nthreads;J++)
+        {
+            *rseed = nextfunc(*rseed);
+            rseeds.data[J] = *rseed + 13*J;
+        }
+        for(J=0;J<nthreads;J++)
+        {
+            threads[J] = new(std::nothrow) std::thread(
+                &amsarray_rand_threadf2<T,T2>,
+                out,
+                &rseeds,
+                randfunc,
+                low,highexcl,
+                J,nthreads
+            );
+            if(threads[J]==NULL)
+            {
+                //handle thread allocation error
+                printf("warning: amsarray_random_threadexec2:: thread %d failed to allocate.\n",J);
+            }
+        }
+        for(J=0;J<nthreads;J++)
+        {
+            if(threads[J]!=NULL)
+            { 
+                threads[J]->join(); 
+                delete threads[J]; 
+                threads[J]= NULL;
+            }
+        }
+    }
+
+    return;
+}
+
+amsarray<int> amsarray_randint(amsarray_size_t N, int low, int highexcl,
+    amsmu_randt1 *rseed)
+{
+    amsarray<int> ret;
+
+    if(N<=0)
+    {
+        ret.resize(0);
+        return ret;
+    }
+
+    amsarray_random_threadexec2(
+        &ret,
+        N,
+        &(randint),
+        &(seed32_next),
+        low,highexcl,
+        rseed
+    );
+
+    return ret;
+}
+
+amsarray<int64_t> amsarray_randintl(amsarray_size_t N, int64_t low, int64_t highexcl, amsmu_randt2 *rseed)
+{
+    amsarray<int64_t> ret;
+
+    if(N<=0)
+    {
+        ret.resize(0);
+        return ret;
+    }
+
+    amsarray_random_threadexec2(
+        &ret,
+        N,
+        &(randintl),
+        &(seed64_next),
+        low,highexcl,
+        rseed
+    );
+
+    return ret;
+}
+
+}; //end namespace rand
+}; //end namespace ams

src/amsmathutil25/testing/amsmathtuil25_test1.cpp

@@ -87,6 +87,104 @@ void test_amsarray2()
     {
         printf("q[%ld]=%1.3f\n",I,q[I]);
     }
+}
+
+void test_amsarray_select()
+{
+    int I;
+    amsarray<double> a;
+    amsarray<amsarray_size_t> b;
+    amsarray<double> c;
+
+    a = ams::rand::amsarray_rand(100);
+
+    for(I=0;I<10;I++)
+    {
+        printf("a[%d]=%1.3f\n",I,a[I]);
+    }
+
+    b = amsarray<amsarray_size_t>({1,3,5,7});
+
+    c = a.select(b);
+
+    for(I=0;I<b.length;I++)
+    {
+        printf("b[%d] = %d: a[%d]=%1.3f, c[%d] = %1.3f\n",I,(int)b[I],(int)b[I],a[b[I]],I,c[I]);
+    }
+
+    a = ams::rand::amsarray_rand(30000);
+    b = (amsarray<amsarray_size_t>)ams::rand::amsarray_randintl(10000,0,a.length);
+    c = a.select(b);
+
+    for(I=b.length/2;I<b.length/2+10;I++)
+    {
+        printf("b[%d] = %d: a[%d]=%1.3f, c[%d] = %1.3f\n",I,(int)b[I],(int)b[I],a[b[I]],I,c[I]);
+    }
+
+    for(I=b.length-10;I<b.length;I++)
+    {
+        printf("b[%d] = %d: a[%d]=%1.3f, c[%d] = %1.3f\n",I,(int)b[I],(int)b[I],a[b[I]],I,c[I]);
+    }
+
+
+
+    return;
+}
+
+void test_amsarray_sort1()
+{
+    int I;
+    amsarray<double> a,b;
+    amsarray<amsarray_size_t> q;
+    int K;
+    //int J;
+    amsarray_size_t bs;
+    double t0,t1;
+
+    ams::rand::seed_init_timer();
+
+
+    for(K=0;K<10;K++)
+    {
+        a = ams::rand::amsarray_rand(10000);
+        //_debug_amsarray_print(&a,1);
+        a.sort();
+        //_debug_amsarray_print(&a,1);
+
+        bool ordertest = 1;
+        for(I=0;I<a.length-1;I++)
+        {
+            if(a[I]>a[I+1])
+            {
+                ordertest = 0;
+                break;
+            }
+        }
+        if(ordertest==1)
+            printf("ordering test: passed.\n");
+        else
+            printf("ordering test: failed.\n");
+    }
+
+    bs = 100;
+    for(K=0;K<10;K++)
+    {
+        b = ams::rand::amsarray_rand(bs);
+        t0 = time_msec();
+        b.sort();
+        t1 = time_msec();
+        printf("sorted %ld in %1.3f msec.\n",bs,t1-t0);
+        bs*=2;
+    };
+    //a.print();
+    
+    // q = permutation_identity(10000);
+
+
+    for(I=0;I<a.length && I<10;I++)
+    {
+        printf("a[%d]=%1.5f\n",I,a[I]);
+    }
 
 
 }

src/amsmathutil25/util/amsmathutil25_amsarray.cpp

@@ -6,7 +6,6 @@ namespace ams
 template<> void amsarray<int>::print(bool newline,int printstyle)
 {
     amsarray_size_t I;
-
     printf("{");
     if(data!=NULL)
     {
@@ -44,7 +43,6 @@ template<> void amsarray<long>::print(bool newline,int printstyle)
 template<> void amsarray<float>::print(bool newline,int printstyle)
 {
     amsarray_size_t I;
-
     printf("{");
     if(data!=NULL)
     {
@@ -63,7 +61,6 @@ template<> void amsarray<float>::print(bool newline,int printstyle)
 template<> void amsarray<double>::print(bool newline,int printstyle)
 {
     amsarray_size_t I;
-
     printf("{");
     if(data!=NULL)
     {
@@ -79,10 +76,49 @@ template<> void amsarray<double>::print(bool newline,int printstyle)
     return;
 }
 
+void _debug_amsarray_print(amsarray<long> *array, bool newline,int printstyle)
+{
+    amsarray_size_t I;
+
+    printf("{");
+    if(array->data!=NULL)
+    {
+        for(I=0;I<array->length-1;I++)
+        {
+            printf("%ld,",array->data[I]);
+        }
+        if(array->length>0) printf("%ld",array->data[array->length-1]);
+    }
+    printf("}");
+    if(newline==1) printf("\n");
+
+    return;
+}
+
+
+void _debug_amsarray_print(amsarray<double> *array, bool newline,int printstyle)
+{
+    amsarray_size_t I;
+
+    printf("{");
+    if(array->data!=NULL)
+    {
+        for(I=0;I<array->length-1;I++)
+        {
+            printf("%1.4f,",array->data[I]);
+        }
+        if(array->length>0) printf("%1.4f",array->data[array->length-1]);
+    }
+    printf("}");
+    if(newline==1) printf("\n");
+
+    return;
+}
+
 // Explicit Class Instantiations?
-template class amsarray<int>; 
-template class amsarray<float>; 
-template class amsarray<int64_t>; 
-template class amsarray<double>; 
+// template class amsarray<int>; 
+// template class amsarray<float>; 
+// template class amsarray<int64_t>; 
+// template class amsarray<double>; 
 
 };

src/amsmathutil25/util/amsmathutil25_amsarray_sort.cpp

@@ -63,8 +63,8 @@ void amsarray_permutation_swap(amsarray<amsarray_size_t> *permutation, amsarray_
 {
     amsarray_size_t tmp;
     tmp = permutation->data[I];
-    permutation->data[J] = permutation->data[I];
-    permutation->data[I] = tmp;
+    permutation->data[I] = permutation->data[J];
+    permutation->data[J] = tmp;
     return;
 }
 

src/amsmathutil25/util/amsmathutil25_utils1.cpp

@@ -0,0 +1,16 @@
+#include <amsmathutil25/amsmathutil25.hpp>
+
+
+namespace ams
+{
+
+//returns time in msec
+double time_msec()
+{
+    double msec = 0.0;
+    clock_t c = clock();
+    msec = (double)c/(double) CLOCKS_PER_SEC * 1000.0;
+    return msec;
+} 
+
+}; //end namespace ams

src/main.cpp

@@ -7,5 +7,8 @@ int main(int argc, char* argv[])
 
     //ams::amsmathutil25::test_amsarray1();
     //ams::amsmathutil25::test_amsarray2();
+    //ams::amsmathutil25::test_amsarray_select();
+    ams::amsmathutil25::test_amsarray_sort1();
+    
     return ret;
 }