updates

build/make.linux64.test.py

@@ -16,7 +16,7 @@ builddir = "./build_linux64"
 doinstall = True   #copies the build_output to the install dir when finished
 cc = "g++" #compiler
 cflags = "-fPIC"
-libraries = "-l{}".format(libname)
+libraries = "-l{} -lamsmathutil2.linux64".format(libname)
 libdirs = "-L{} -L{}/lib".format(builddir,commondir)
 linkerflags = "-static -static-libgcc -Wl,-rpath=."
 srcexts = [".c",".cpp"]

build/make.mingw64.test.py

@@ -17,7 +17,7 @@ builddir = "./build_mingw64"
 doinstall = False   #copies the build_output to the install dir when finished
 cc = "x86_64-w64-mingw32-g++" #compiler
 cflags = "-fPIC -O3"
-libraries = "-l{}".format(libname)
+libraries = "-l{} -lamsmathutil2.mingw64".format(libname)
 libdirs = "-L{} -L{}/lib -L{}/lib".format(builddir,commondir,depdir)
 linkerflags = "-static -static-libgcc -Wl,-rpath=."
 srcexts = [".c",".cpp"]

include/amscppnarray/amscppnarray.hpp

@@ -4,11 +4,16 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include <time.h>
+
 #include <vector>
 #include <new>
 #include <initializer_list>
 #include <thread>
 
+#include <amsmath/amsmath.hpp>
+#include <amsutil/amsutil.hpp>
+
 namespace ams
 {
 namespace narray
@@ -114,5 +119,6 @@ void test_narray3();
 #include <amscppnarray/amscppnarray_implinits.hpp>
 
 #include <amscppnarray/amscppnarray_random.hpp>
+#include <amscppnarray/amscppnarray_math.hpp>
 
 #endif

include/amscppnarray/amscppnarray_math.hpp

@@ -0,0 +1,28 @@
+#ifndef __AMSCPPNARRAY_MATH_HPP__
+#define __AMSCPPNARRAY_MATH_HPP__
+
+namespace ams
+{
+namespace narray
+{
+
+template<typename T> T narray_max(const narray<T> &q);
+template<typename T> T narray_min(const narray<T> &q);
+template<typename T> double narray_mean(const narray<T> &q);
+template<typename T> double narray_var(const narray<T> &q);
+template<typename T> double narray_stdev(const narray<T> &q);
+template<typename T> double narray_L2norm(const narray<T> &q);
+
+template<typename T> T narray_sum(const narray<T> &q);
+
+template<typename T> narray<T> narray_abs(narray<T> &q);
+
+void narray_testmath1();
+void narray_testmath2();
+
+}; //end namespaces
+};
+
+#include <amscppnarray/amscppnarray_mathimpl.hpp>
+
+#endif

include/amscppnarray/amscppnarray_mathimpl.hpp

@@ -0,0 +1,381 @@
+#ifndef __AMSCPPNARRAY_MATHIMPL_HPP__
+#define __AMSCPPNARRAY_MATHIMPL_HPP__
+
+namespace ams
+{
+namespace narray
+{
+
+    template<typename T> T narray_max(const narray<T> &q)
+    {
+        return q.max();
+    }
+
+    template<typename T> T narray_min(const narray<T> &q)
+    {
+        return q.min();
+    }
+
+    template<typename T> void narray_sum_tf(
+        const narray<T> *in,
+        narray<T> *sums,
+        int threadnum,
+        int nthreads
+    )
+    {
+        int I0,I1,Is,I;
+    
+        Is = in->length/nthreads;
+        I0 = Is*threadnum;
+        I1 = (threadnum>=(nthreads-1))? in->length: Is*(threadnum+1);
+        for(I=I0;I<I1;I++)
+        {
+            sums->data[threadnum] = sums->data[threadnum] + in->data[I];
+        }
+    }
+
+    template<typename T> T narray_sum(const narray<T> &q)
+    {
+        T ret = T();
+        narray<T> sums;
+        narray_size_t I;
+        int J;
+        int nthreads;
+        std::vector<std::thread*> threads;
+
+        if(q.length<=0)
+        {
+            return ret;
+        }
+
+        if(q.length<narray_thread_sz)
+        {
+            ret = T();
+            for(I=0;I<q.length;I++)
+            {
+                ret = ret + q.data[I];
+            }
+        }
+        else
+        {
+            //threaded operation
+            
+            nthreads = std::thread::hardware_concurrency();
+            nthreads = (nthreads<1) ? 1 : nthreads;
+            nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+            threads.resize(nthreads);
+            sums.resize(nthreads);
+            sums.setall(T());
+            for(J=0;J<nthreads;J++)
+            {
+                threads[J] = new(std::nothrow) std::thread(
+                    &narray_sum_tf<T>,
+                    &q,
+                    &sums,
+                    J,nthreads
+                );
+                if(threads[J]==NULL)
+                {
+                    printf("narray_sum: warning: thread %d failed to initialize.\n",J);
+                    //handle errors
+                }
+            }
+            for(J=0;J<nthreads;J++)
+            {
+                if(threads[J]!=NULL)
+                { 
+                    threads[J]->join(); 
+                    delete threads[J]; 
+                    threads[J]= NULL;
+                }
+            }
+            ret = T();
+            for(J=0;J<nthreads;J++)
+            {
+                ret = ret + sums.data[J];
+            }
+        }
+        return ret;
+    }
+
+    template<typename T> double narray_mean(const narray<T> &q)
+    {
+        double ret = 0.0;
+        T sum = narray_sum(q);
+        ret = (double)sum/((double)q.length);
+        return ret;
+    }
+
+    template<typename T> void narray_variance_tf(
+        const narray<T> *in,
+        narray<double> *sums,
+        double mean,
+        int threadnum,
+        int nthreads
+    )
+    {
+        int I0,I1,Is,I;
+        double q2;
+        Is = in->length/nthreads;
+        I0 = Is*threadnum;
+        I1 = (threadnum>=(nthreads-1))? in->length: Is*(threadnum+1);
+        for(I=I0;I<I1;I++)
+        {
+            q2 = in->data[I]-mean;
+            sums->data[threadnum] = sums->data[threadnum] + q2*q2/((double)in->length);
+        }
+    }
+
+    template<typename T> double narray_var(const narray<T> &q)
+    {
+        double ret = 0.0;
+        narray<double> sums;
+        narray_size_t I;
+        int J;
+        int nthreads;
+        std::vector<std::thread*> threads;
+
+        double mn = ams::narray::narray_mean(q);
+        double q2;
+
+        if(q.length<=0)
+        {
+            return ret;
+        }
+
+        if(q.length<narray_thread_sz)
+        {
+            ret = 0.0;
+            for(I=0;I<q.length;I++)
+            {
+                q2 = (double)q.data[I]-mn;
+                ret = ret + q2*q2/((double)q.length);
+            }
+        }
+        else
+        {
+            //threaded operation
+            
+            nthreads = std::thread::hardware_concurrency();
+            nthreads = (nthreads<1) ? 1 : nthreads;
+            nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+            threads.resize(nthreads);
+            sums.resize(nthreads);
+            sums.setall(T());
+            for(J=0;J<nthreads;J++)
+            {
+                threads[J] = new(std::nothrow) std::thread(
+                    &narray_variance_tf<T>,
+                    &q,
+                    &sums,
+                    mn,
+                    J,nthreads
+                );
+                if(threads[J]==NULL)
+                {
+                    printf("narray_var: warning: thread %d failed to initialize.\n",J);
+                    //handle errors
+                }
+            }
+            for(J=0;J<nthreads;J++)
+            {
+                if(threads[J]!=NULL)
+                { 
+                    threads[J]->join(); 
+                    delete threads[J]; 
+                    threads[J]= NULL;
+                }
+            }
+            ret = T();
+            for(J=0;J<nthreads;J++)
+            {
+                ret = ret + sums.data[J];
+            }
+        }
+        return ret;
+    }
+
+    template<typename T> double narray_stdev(const narray<T> &q)
+    {
+        double ret = narray_var(q);
+        ret = ::sqrt(ret);
+        return ret;
+    }
+
+    template<typename T> void narray_L2norm_tf(
+        const narray<T> *in,
+        narray<double> *sums,
+        int threadnum,
+        int nthreads
+    )
+    {
+        int I0,I1,Is,I;
+        double q2;
+        Is = in->length/nthreads;
+        I0 = Is*threadnum;
+        I1 = (threadnum>=(nthreads-1))? in->length: Is*(threadnum+1);
+        for(I=I0;I<I1;I++)
+        {
+            q2 = in->data[I];
+            sums->data[threadnum] = sums->data[threadnum] + q2*q2;
+        }
+    }
+
+    template<typename T> double narray_L2norm(const narray<T> &q)
+    {
+        double ret = 0.0;
+        narray<double> sums;
+        narray_size_t I;
+        int J;
+        int nthreads;
+        std::vector<std::thread*> threads;
+        double q2;
+
+        if(q.length<=0)
+        {
+            return ret;
+        }
+
+        if(q.length<narray_thread_sz)
+        {
+            ret = 0.0;
+            for(I=0;I<q.length;I++)
+            {
+                q2 = (double)q.data[I];
+                ret = ret + q2*q2;
+            }
+        }
+        else
+        {
+            //threaded operation
+            
+            nthreads = std::thread::hardware_concurrency();
+            nthreads = (nthreads<1) ? 1 : nthreads;
+            nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+            threads.resize(nthreads);
+            sums.resize(nthreads);
+            sums.setall(T());
+            for(J=0;J<nthreads;J++)
+            {
+                threads[J] = new(std::nothrow) std::thread(
+                    &narray_L2norm_tf<T>,
+                    &q,
+                    &sums,
+                    J,nthreads
+                );
+                if(threads[J]==NULL)
+                {
+                    printf("narray_L2norm: warning: thread %d failed to initialize.\n",J);
+                    //handle errors
+                }
+            }
+            for(J=0;J<nthreads;J++)
+            {
+                if(threads[J]!=NULL)
+                { 
+                    threads[J]->join(); 
+                    delete threads[J]; 
+                    threads[J]= NULL;
+                }
+            }
+            ret = T();
+            for(J=0;J<nthreads;J++)
+            {
+                ret = ret + sums.data[J];
+            }
+        }
+
+        ret = ::sqrt(ret);
+        return ret;
+    }
+
+
+    template<typename T> void narray_abs_tf(
+        const narray<T> *in,
+        narray<T> *out,
+        int threadnum,
+        int nthreads
+    )
+    {
+        int I0,I1,Is,I;
+        double q2;
+        Is = in->length/nthreads;
+        I0 = Is*threadnum;
+        I1 = (threadnum>=(nthreads-1))? in->length: Is*(threadnum+1);
+        for(I=I0;I<I1;I++)
+        {
+            out->data[I] = (in->data[I]<T(0)) ? -in->data[I] : in->data[I];
+        }
+    }
+
+
+    template<typename T> narray<T> narray_abs(narray<T> &q)
+    {
+        narray<T> ret;
+        narray_size_t I;
+        int J;
+        int nthreads;
+        std::vector<std::thread*> threads;
+
+        if(q.length<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        ret.resize(q.length);
+        if(ret.length!=q.length)
+        {
+            printf("narray_abs: error - could not allocate return array.\n");
+            ret.resize(0);
+            return ret;
+        }
+
+        if(q.length<narray_thread_sz)
+        {
+            for(I=0;I<q.length;I++)
+            {
+                ret.data[I] = (q.data[I]<T(0)) ? -q.data[I] : q.data[I];
+            }
+        }
+        else
+        {
+            //threaded operation
+            
+            nthreads = std::thread::hardware_concurrency();
+            nthreads = (nthreads<1) ? 1 : nthreads;
+            nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+            threads.resize(nthreads);
+            for(J=0;J<nthreads;J++)
+            {
+                threads[J] = new(std::nothrow) std::thread(
+                    &narray_abs_tf<T>,
+                    &q,
+                    &ret,
+                    J,nthreads
+                );
+                if(threads[J]==NULL)
+                {
+                    printf("narray_abs: warning: thread %d failed to initialize.\n",J);
+                    //handle errors
+                }
+            }
+            for(J=0;J<nthreads;J++)
+            {
+                if(threads[J]!=NULL)
+                { 
+                    threads[J]->join(); 
+                    delete threads[J]; 
+                    threads[J]= NULL;
+                }
+            }
+        }
+
+        return ret;
+    }
+
+    
+}; //end namespaces
+};
+
+#endif

include/amscppnarray/amscppnarray_random.hpp

@@ -4,8 +4,71 @@
 namespace ams
 {
 namespace narray
+{
+namespace rand
 {
 
+    //Random functions for numeric array class
+
+    ////////////////////////////////////////////////////////////////////////
+    //Simple Deterministic Psuedorandom Number Generator (32 bit version) //
+    ////////////////////////////////////////////////////////////////////////
+
+    static const int32_t dpr32_mod = ( ((int32_t)1) << ((int32_t)30) ) - (int32_t)1;
+    static const int32_t dpr32_mult1 = ( (int32_t) 1201633 );
+    static int32_t dpr32_rseed = 0; // global random number seed for numeric array class
+    
+    int32_t dpr32_nextseed(int32_t seed);
+    double dpr32_randd(int32_t *seed);
+    float dpr32_randf(int32_t *seed);
+    double dpr32_gaussian(int32_t *seed);
+    float dpr32_gaussianf(int32_t *seed);
+    int dpr32_randint(int32_t *seed, int low, int high);
+    
+    ams::vect2 dpr32_randuvect2(int32_t *seed);
+    ams::vect3 dpr32_randuvect3(int32_t *seed);
+    ams::vect4 dpr32_randuvect4(int32_t *seed);
+    
+    ams::vect2f dpr32_randuvect2f(int32_t *seed);
+    ams::vect3f dpr32_randuvect3f(int32_t *seed);
+    ams::vect4f dpr32_randuvect4f(int32_t *seed);
+
+    ////////////////////////////////////////////
+    //Numeric array threaded random generators//
+    ////////////////////////////////////////////
+
+    int32_t get_rseed();
+    void set_rseed(int32_t seed);
+    void set_rseed_withtimer();
+
+    narray<double> narray_rand(narray_size_t N, 
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    narray<float> narray_randf(narray_size_t N,
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    narray<double> narray_randgauss(narray_size_t N,
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    narray<float> narray_randgaussf(narray_size_t N,
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    
+    narray<int> narray_randint(narray_size_t N, int low, int highexcl,
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    narray<vect2> narray_randuvect2(narray_size_t N,
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    narray<vect3> narray_randuvect3(narray_size_t N,
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    narray<vect4> narray_randuvect4(narray_size_t N,
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    narray<vect2f> narray_randuvect2f(narray_size_t N,
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    narray<vect3f> narray_randuvect3f(narray_size_t N,
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    narray<vect4f> narray_randuvect4f(narray_size_t N,
+        int32_t *rseed = &(ams::narray::rand::dpr32_rseed));
+    
+
+    void test_random1();
+    
+};
 };
 };
 

src/amscppnarray/amscppnarray_math.cpp

@@ -0,0 +1,85 @@
+#include <amscppnarray/amscppnarray.hpp>
+
+namespace ams
+{
+namespace narray
+{
+
+    static void _intl_printarr(narray<int> &q)
+    {
+        int I;
+        if(q.length>0)
+        {
+            printf("{");
+            for(I=0;I<q.length-1;I++) printf("%d,",q[I]);
+            printf("%d}",q[q.length-1]);
+        }
+        else
+        {
+            printf("{}");
+        }
+    }
+
+    static void _intl_printarr(narray<double> &q)
+    {
+        int I;
+        if(q.length>0)
+        {
+            printf("{");
+            for(I=0;I<q.length-1;I++) printf("%1.3f,",q[I]);
+            printf("%1.3f}",q[q.length-1]);
+        }
+        else
+        {
+            printf("{}");
+        }
+    }
+
+    void narray_testmath1()
+    {
+        narray<int> a = {-4,-3,-2,-1,0,1,2,3,4,5};
+        narray<int> b;
+        int asum = narray_sum(a);
+        double amean = narray_mean(a);
+
+        printf("sum(a)=%d\n",asum);
+        printf("mean(a)=%1.3f\n",amean);
+        printf("stdev(a)=%1.3f\n",narray_stdev(a));
+        printf("L2norm(a)=%1.3f\n",narray_L2norm(a));
+        printf("min(a)=%d\n",narray_min(a));
+        printf("max(a)=%d\n",narray_max(a));
+
+        b = narray_abs(a);
+        _intl_printarr(b); printf("\n");
+
+        return;
+    }
+
+    void narray_testmath2()
+    {
+        rand::set_rseed_withtimer();
+
+        narray<double> unif = ams::narray::rand::narray_rand(10000);
+        double mn1 = narray_mean(unif);
+        double std1 = narray_stdev(unif);
+        double mn2,std2;
+        
+        
+
+        //_intl_printarr(unif); printf("\n");
+
+        unif = unif*2.0-1.0;
+        //_intl_printarr(unif); printf("\n");
+        mn2 = narray_mean(unif);
+        std2 = narray_stdev(unif);
+
+
+        printf("mn1=%1.3f std1=%1.3f\n",mn1,std1);
+        printf("mn2=%1.3f std2=%1.3f\n",mn2,std2);
+        
+
+        return;
+    }
+
+};
+};

src/amscppnarray/amscppnarray_random.cpp

@@ -3,9 +3,584 @@
 namespace ams
 {
 namespace narray
+{
+namespace rand
 {
 
-
+    int32_t dpr32_nextseed(int32_t seed)
+    {
+        int32_t sret = seed;
+        sret = ams::mod(sret*dpr32_mult1+1,dpr32_mod);
+        return sret;
+    }
     
+    double dpr32_randd(int32_t *seed)
+    {
+        double ret;
+        *seed = ams::narray::rand::dpr32_nextseed(*seed);
+        ret = (double)*seed/(double)(dpr32_mod-1);
+        return ret;
+    }
+    
+    float dpr32_randf(int32_t *seed)
+    {
+        float ret;
+        *seed = ams::narray::rand::dpr32_nextseed(*seed);
+        ret = (float)*seed/(float)(dpr32_mod-1);
+        return ret;
+    }
+    
+    double dpr32_gaussian(int32_t *seed)
+    {
+        double ret = 0.0;
+        double u1,u2;
+        u1 = ams::narray::rand::dpr32_randd(seed);
+        u2 = ams::narray::rand::dpr32_randd(seed);
+        
+        if(u1>0.0)
+        {
+            ret = ::sqrt(-2.0*::log(u1))*::cos(2.0*pi*u2);
+        }
+    
+        return ret;
+    }
+    
+    
+    float dpr32_gaussianf(int32_t *seed)
+    {
+        float ret = 0.0f;
+        float u1,u2;
+        u1 = ams::narray::rand::dpr32_randf(seed);
+        u2 = ams::narray::rand::dpr32_randf(seed);
+        
+        if(u1>0.0)
+        {
+            ret = ::sqrt(-2.0f*::log(u1))*::cos(2.0f*pif*u2);
+        }
+    
+        return ret;
+    }
+    
+    vect2 dpr32_randuvect2(int32_t *seed)
+    {
+        vect2 ret = vect2(0,0);
+        double th = ams::narray::rand::dpr32_randd(seed)*2.0*pi;
+        ret = vect2(cos(th),sin(th));
+        return ret;
+    }
+    
+    vect3 dpr32_randuvect3(int32_t *seed)
+    {
+        vect3 ret;
+        double az,el;
+        el = ::acos(2.0*dpr32_randd(seed)-1.0)-pi/2.0;
+        az = ams::narray::rand::dpr32_randd(seed)*2.0*pi;
+    
+        ret = vect3(cos(az)*cos(el),sin(az)*cos(el),sin(el));
+    
+        return ret;
+    }
+    
+    vect4 dpr32_randuvect4(int32_t *seed)
+    {
+        vect4 ret;
+        double x,y,z,w;
+    
+        x = ams::narray::rand::dpr32_gaussian(seed);
+        y = ams::narray::rand::dpr32_gaussian(seed);
+        z = ams::narray::rand::dpr32_gaussian(seed);
+        w = ams::narray::rand::dpr32_gaussian(seed);
+        ret = vect4(x,y,z,w);
+        ret = vnormalize(ret);
+    
+        return ret;
+    }
+    
+    vect2f dpr32_randuvect2f(int32_t *seed)
+    {
+        vect2f ret = vect2f(0,0);
+        float th = ams::narray::rand::dpr32_randf(seed)*2.0*pi;
+        ret = vect2f(cos(th),sin(th));
+        return ret;
+    }
+    
+    vect3f dpr32_randuvect3f(int32_t *seed)
+    {
+        vect3f ret;
+        float az,el;
+        el = ::acos(2.0f*dpr32_randf(seed)-1.0)-pif/2.0;
+        az = ams::narray::rand::dpr32_randf(seed)*2.0f*pif;
+    
+        ret = vect3f(cos(az)*cos(el),sin(az)*cos(el),sin(el));
+    
+        return ret;
+    }
+    
+    vect4f dpr32_randuvect4f(int32_t *seed)
+    {
+        vect4f ret;
+        float x,y,z,w;
+    
+        x = ams::narray::rand::dpr32_gaussianf(seed);
+        y = ams::narray::rand::dpr32_gaussianf(seed);
+        z = ams::narray::rand::dpr32_gaussianf(seed);
+        w = ams::narray::rand::dpr32_gaussianf(seed);
+        ret = vect4f(x,y,z,w);
+        ret = vnormalize(ret);
+    
+        return ret;
+    }
+
+    int32_t get_rseed()
+    {
+        return ams::narray::rand::dpr32_rseed;
+    }
+
+    void set_rseed(int32_t seed)
+    {
+        ams::narray::rand::dpr32_rseed = seed;
+    }
+
+    void set_rseed_withtimer()
+    {
+        int32_t t1 = (int32_t)time(NULL);
+        int32_t t2 = (int32_t)(::fmod((double)clock()/((double)CLOCKS_PER_SEC)*1000.0f,36000.0f));
+        ams::narray::rand::dpr32_rseed = (int32_t)t1 + (int32_t)t2;
+    }
+
+    template<typename T> void narray_rand_threadf1(
+        narray<T> *out,
+        narray<int32_t> *rseeds,
+        T (*randfunc)(int32_t *),
+        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 narray_random_threadexec1(
+        narray<T> *out,
+        narray_size_t N,
+        T (*randfunc)(int32_t *),
+        int32_t *rseed
+    )
+    {
+        narray_size_t I;
+        int J;
+        int nthreads;
+        std::vector<std::thread*> threads;
+        narray<int32_t> rseeds;
+        int res;
+
+        res = out->resize(N);
+        if(res!=narray_success)
+        {
+            out->resize(0);
+            return;
+        }
+
+        if(N<narray_thread_sz)
+        {
+            //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>narray_max_threads) ? narray_max_threads : nthreads;
+            threads.resize(nthreads);
+            rseeds.resize(nthreads);
+            for(J=0;J<nthreads;J++)
+            {
+                *rseed = dpr32_nextseed(*rseed);
+                rseeds.data[J] = *rseed + 13*J;
+            }
+            for(J=0;J<nthreads;J++)
+            {
+                threads[J] = new(std::nothrow) std::thread(
+                    &narray_rand_threadf1<T>,
+                    out,
+                    &rseeds,
+                    randfunc,
+                    J,nthreads
+                );
+                if(threads[J]==NULL)
+                {
+                    //handle thread allocation error
+                    printf("warning: narray_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;
+    }
+
+    narray<double> narray_rand(narray_size_t N,
+        int32_t *rseed)
+    {
+        narray<double> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec1(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_randd),
+            rseed
+        );
+
+        return ret;
+    }
+
+    narray<float> narray_randf(narray_size_t N,
+        int32_t *rseed)
+    {
+        narray<float> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec1(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_randf),
+            rseed
+        );
+
+        return ret;
+    }
+
+    narray<double> narray_randgauss(narray_size_t N,
+        int32_t *rseed)
+    {
+        narray<double> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec1(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_gaussian),
+            rseed
+        );
+
+        return ret;
+    }
+
+    narray<float> narray_randgaussf(narray_size_t N,
+        int32_t *rseed)
+    {
+        narray<float> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec1(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_gaussianf),
+            rseed
+        );
+
+        return ret;
+    }
+
+    narray<vect2> narray_randuvect2(narray_size_t N,
+        int32_t *rseed)
+    {
+        narray<vect2> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec1(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_randuvect2),
+            rseed
+        );
+
+        return ret;
+    }
+
+    narray<vect3> narray_randuvect3(narray_size_t N,
+        int32_t *rseed)
+    {
+        narray<vect3> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec1(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_randuvect3),
+            rseed
+        );
+
+        return ret;
+    }
+
+    narray<vect4> narray_randuvect4(narray_size_t N,
+        int32_t *rseed)
+    {
+        narray<vect4> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec1(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_randuvect4),
+            rseed
+        );
+
+        return ret;
+    }
+
+    narray<vect2f> narray_randuvect2f(narray_size_t N,
+        int32_t *rseed)
+    {
+        narray<vect2f> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec1(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_randuvect2f),
+            rseed
+        );
+
+        return ret;
+    }
+
+    narray<vect3f> narray_randuvect3f(narray_size_t N,
+        int32_t *rseed)
+    {
+        narray<vect3f> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec1(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_randuvect3f),
+            rseed
+        );
+
+        return ret;
+    }
+
+    narray<vect4f> narray_randuvect4f(narray_size_t N,
+        int32_t *rseed)
+    {
+        narray<vect4f> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec1(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_randuvect4f),
+            rseed
+        );
+
+        return ret;
+    }
+
+
+    int dpr32_randint(int32_t *seed, int low, int high)
+    {
+        int ret = 0;
+        int val;
+
+        if(high-low>0)
+        {
+            *seed = dpr32_nextseed(*seed);
+            val = low + (int)(*seed)%(high-low);
+        }
+
+        return ret;
+    }
+
+    template<typename T> void narray_rand_threadf2(
+        narray<T> *out,
+        narray<int32_t> *rseeds,
+        T (*randfunc)(int32_t *, T, T),
+        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(&(rseeds->data[threadnum]),low,highexcl);
+        }
+        return;
+    }
+
+    template<typename T> void narray_random_threadexec2(
+        narray<T> *out,
+        narray_size_t N,
+        T (*randfunc)(int32_t *, T, T),
+        T low,
+        T highexcl,
+        int32_t *rseed
+    )
+    {
+        narray_size_t I;
+        int J;
+        int nthreads;
+        std::vector<std::thread*> threads;
+        narray<int32_t> rseeds;
+        int res;
+
+        res = out->resize(N);
+        if(res!=narray_success)
+        {
+            out->resize(0);
+            return;
+        }
+
+        if(N<narray_thread_sz)
+        {
+            //single threaded
+            for(I=0;I<N;I++)
+            {
+                out->data[I] = randfunc(rseed,low,highexcl);
+            }
+        }
+        else
+        {
+            //multi-threaded operation
+            nthreads = std::thread::hardware_concurrency();
+            nthreads = (nthreads<1) ? 1 : nthreads;
+            nthreads = (nthreads>narray_max_threads) ? narray_max_threads : nthreads;
+            threads.resize(nthreads);
+            rseeds.resize(nthreads);
+            for(J=0;J<nthreads;J++)
+            {
+                *rseed = dpr32_nextseed(*rseed);
+                rseeds.data[J] = *rseed + 13*J;
+            }
+            for(J=0;J<nthreads;J++)
+            {
+                threads[J] = new(std::nothrow) std::thread(
+                    &narray_rand_threadf2<T>,
+                    out,
+                    &rseeds,
+                    randfunc,
+                    low,highexcl,
+                    J,nthreads
+                );
+                if(threads[J]==NULL)
+                {
+                    //handle thread allocation error
+                    printf("warning: narray_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;
+    }
+
+    narray<int> narray_randint(narray_size_t N, int low, int highexcl,
+        int32_t *rseed)
+    {
+        narray<int> ret;
+
+        if(N<=0)
+        {
+            ret.resize(0);
+            return ret;
+        }
+
+        narray_random_threadexec2(
+            &ret,
+            N,
+            &(ams::narray::rand::dpr32_randint),
+            low,highexcl,
+            rseed
+        );
+
+        return ret;
+    }
+
+
+}; //end namespace rand
 }; //end namespaces
 };

src/main.cpp

@@ -7,6 +7,7 @@ int main(int argc, char* argv[])
     //ams::narray::test_narray1();
     //ams::narray::test_narray2();
     //ams::narray::test_narray3();
+    ams::narray::narray_testmath2();
 
     return ret;
 }