diff --git a/build_linux64/libamsculib3.linux64.a b/build_linux64/libamsculib3.linux64.a
index e06ae52..3905e0b 100644
Binary files a/build_linux64/libamsculib3.linux64.a and b/build_linux64/libamsculib3.linux64.a differ
diff --git a/build_linux64/objstore/amscu_comp128.o b/build_linux64/objstore/amscu_comp128.o
index 2c90406..0e0fd51 100644
Binary files a/build_linux64/objstore/amscu_comp128.o and b/build_linux64/objstore/amscu_comp128.o differ
diff --git a/build_linux64/objstore/amscu_comp64.o b/build_linux64/objstore/amscu_comp64.o
index ba9bed6..252e8da 100644
Binary files a/build_linux64/objstore/amscu_comp64.o and b/build_linux64/objstore/amscu_comp64.o differ
diff --git a/build_linux64/objstore/amscu_cudafunctions.o b/build_linux64/objstore/amscu_cudafunctions.o
index 1481357..cc8e081 100644
Binary files a/build_linux64/objstore/amscu_cudafunctions.o and b/build_linux64/objstore/amscu_cudafunctions.o differ
diff --git a/build_linux64/objstore/amscu_random.o b/build_linux64/objstore/amscu_random.o
index 17da646..2ba0e79 100644
Binary files a/build_linux64/objstore/amscu_random.o and b/build_linux64/objstore/amscu_random.o differ
diff --git a/build_linux64/objstore/amscuarray.o b/build_linux64/objstore/amscuarray.o
index 5ab9e32..df0bbb4 100644
Binary files a/build_linux64/objstore/amscuarray.o and b/build_linux64/objstore/amscuarray.o differ
diff --git a/build_linux64/objstore/amscuarray_dops.o b/build_linux64/objstore/amscuarray_dops.o
index 77ceb4b..2ac9ce6 100644
Binary files a/build_linux64/objstore/amscuarray_dops.o and b/build_linux64/objstore/amscuarray_dops.o differ
diff --git a/build_linux64/objstore/amscugeom.o b/build_linux64/objstore/amscugeom.o
index dfe3b7a..2aa1ff0 100644
Binary files a/build_linux64/objstore/amscugeom.o and b/build_linux64/objstore/amscugeom.o differ
diff --git a/build_linux64/objstore/amsculib3.o b/build_linux64/objstore/amsculib3.o
index 4820382..fc6935a 100644
Binary files a/build_linux64/objstore/amsculib3.o and b/build_linux64/objstore/amsculib3.o differ
diff --git a/build_linux64/objstore/amscumath.o b/build_linux64/objstore/amscumath.o
index 4087a4b..2034c0a 100644
Binary files a/build_linux64/objstore/amscumath.o and b/build_linux64/objstore/amscumath.o differ
diff --git a/build_linux64/objstore/amscurarray.o b/build_linux64/objstore/amscurarray.o
index 8f08c5f..2354199 100644
Binary files a/build_linux64/objstore/amscurarray.o and b/build_linux64/objstore/amscurarray.o differ
diff --git a/build_linux64/objstore/cuvec2f.o b/build_linux64/objstore/cuvec2f.o
index e54b9f3..82d3901 100644
Binary files a/build_linux64/objstore/cuvec2f.o and b/build_linux64/objstore/cuvec2f.o differ
diff --git a/build_linux64/objstore/cuvec3f.o b/build_linux64/objstore/cuvec3f.o
index 0f10658..1174501 100644
Binary files a/build_linux64/objstore/cuvec3f.o and b/build_linux64/objstore/cuvec3f.o differ
diff --git a/build_linux64/objstore/cuvec4f.o b/build_linux64/objstore/cuvec4f.o
index 8f97486..6483ef8 100644
Binary files a/build_linux64/objstore/cuvec4f.o and b/build_linux64/objstore/cuvec4f.o differ
diff --git a/build_linux64/test b/build_linux64/test
index a01602f..131155b 100644
Binary files a/build_linux64/test and b/build_linux64/test differ
diff --git a/include/amsculib3/math/amscu_comp128.hpp b/include/amsculib3/math/amscu_comp128.hpp
index a100d1c..df0b73b 100644
--- a/include/amsculib3/math/amscu_comp128.hpp
+++ b/include/amsculib3/math/amscu_comp128.hpp
@@ -18,23 +18,26 @@ namespace cmp
         __host__ __device__ explicit cucomp128(const double &other);
         __host__ __device__ explicit cucomp128(const double &_real, const double &_imag);
 
-        __host__ __device__ cucomp128& operator=(cucomp128& other);
-        __host__ __device__ const cucomp128& operator=(const cucomp128& other);
-        __host__ __device__ cucomp128& operator=(double& other);
-        __host__ __device__ const cucomp128& operator=(const double& other);
+        __host__ __device__ cucomp128& operator=(const cucomp128& other);
+        __host__ __device__ cucomp128& operator=(const double& other);
         
         __host__ __device__ double& operator[](int& ind);
         __host__ __device__ const double& operator[](const int& ind) const;
 
-        __host__ __device__ cucomp128 operator+(const cucomp128& z);
-        __host__ __device__ cucomp128 operator-(const cucomp128& z);
-        __host__ __device__ cucomp128 operator*(const cucomp128& z);
-        __host__ __device__ cucomp128 operator/(const cucomp128& z);
+        __host__ __device__ cucomp128 operator+(const cucomp128& z) const;
+        __host__ __device__ cucomp128 operator-(const cucomp128& z) const;
+        __host__ __device__ cucomp128 operator*(const cucomp128& z) const;
+        __host__ __device__ cucomp128 operator/(const cucomp128& z) const;
 
-        __host__ __device__ cucomp128 operator+(const double& z);
-        __host__ __device__ cucomp128 operator-(const double& z);
-        __host__ __device__ cucomp128 operator*(const double& z);
-        __host__ __device__ cucomp128 operator/(const double& z);
+        __host__ __device__ friend cucomp128 operator+(const cucomp128& z1, const double& z2);
+        __host__ __device__ friend cucomp128 operator-(const cucomp128& z1, const double& z2);
+        __host__ __device__ friend cucomp128 operator*(const cucomp128& z1, const double& z2);
+        __host__ __device__ friend cucomp128 operator/(const cucomp128& z1, const double& z2);
+
+        __host__ __device__ friend cucomp128 operator+(const double& z1, const cucomp128& z2);
+        __host__ __device__ friend cucomp128 operator-(const double& z1, const cucomp128& z2);
+        __host__ __device__ friend cucomp128 operator*(const double& z1, const cucomp128& z2);
+        __host__ __device__ friend cucomp128 operator/(const double& z1, const cucomp128& z2);
 
         __host__ __device__ friend cucomp128 operator-(const cucomp128& z); //negation sign
 
@@ -54,30 +57,42 @@ namespace cmp
         __host__ __device__ bool iszero() const;
         __host__ __device__ double arg() const;
         __host__ __device__ double mag() const;
+        __host__ __device__ double abs() const;
         __host__ __device__ cucomp128 conj() const;
+
+        //accumulation operators
+        __host__ __device__ cucomp128& operator+=(const cucomp128& z);
+        __host__ __device__ cucomp128& operator-=(const cucomp128& z);
+        __host__ __device__ cucomp128& operator*=(const cucomp128& z);
+        __host__ __device__ cucomp128& operator/=(const cucomp128& z);
+
+        __host__ __device__ cucomp128& operator+=(const double& z);
+        __host__ __device__ cucomp128& operator-=(const double& z);
+        __host__ __device__ cucomp128& operator*=(const double& z);
+        __host__ __device__ cucomp128& operator/=(const double& z);
+        
     };
 
-    __host__ __device__ double arg(cucomp128 z);
-
-    __host__ __device__ double real(cucomp128 z);
-    __host__ __device__ double imag(cucomp128 z);
-    __host__ __device__ cucomp128 sin(cucomp128 z);
-    __host__ __device__ cucomp128 cos(cucomp128 z);
-    __host__ __device__ cucomp128 tan(cucomp128 z);
-    __host__ __device__ cucomp128 exp(cucomp128 z);
-    __host__ __device__ cucomp128 log(cucomp128 z);
-    __host__ __device__ double abs(cucomp128 z);
-    __host__ __device__ cucomp128 conj(cucomp128 z);
+    __host__ __device__ double arg(const cucomp128 &z);
+    __host__ __device__ double real(const cucomp128 &z);
+    __host__ __device__ double imag(const cucomp128 &z);
+    __host__ __device__ cucomp128 sin(const cucomp128 &z);
+    __host__ __device__ cucomp128 cos(const cucomp128 &z);
+    __host__ __device__ cucomp128 tan(const cucomp128 &z);
+    __host__ __device__ cucomp128 exp(const cucomp128 &z);
+    __host__ __device__ cucomp128 log(const cucomp128 &z);
+    __host__ __device__ double abs(const cucomp128 &z);
+    __host__ __device__ cucomp128 conj(const cucomp128 &z);
 
     // //need hyperbolic trig Functions
-    __host__ __device__ cucomp128 cosh(cucomp128 z);
-    __host__ __device__ cucomp128 sinh(cucomp128 z);
-    __host__ __device__ cucomp128 tanh(cucomp128 z);
+    __host__ __device__ cucomp128 cosh(const cucomp128 &z);
+    __host__ __device__ cucomp128 sinh(const cucomp128 &z);
+    __host__ __device__ cucomp128 tanh(const cucomp128 &z);
 
-    __host__ __device__ cucomp128 pow(cucomp128 z1, cucomp128 z2);
+    __host__ __device__ cucomp128 pow(const cucomp128 &z1, const cucomp128 &z2);
 
     // //returns "complex sign" of complex number - 0, or a unit number with same argument
-    __host__ __device__ cucomp128 csgn(cucomp128 z);
+    __host__ __device__ cucomp128 csgn(const cucomp128 &z);
 
 void test_cucomp128_1();
 
diff --git a/src/amsculib3/math/amscu_comp128.cu b/src/amsculib3/math/amscu_comp128.cu
index 6e1aef8..9645b90 100644
--- a/src/amsculib3/math/amscu_comp128.cu
+++ b/src/amsculib3/math/amscu_comp128.cu
@@ -5,438 +5,530 @@ namespace amscuda
 namespace cmp
 {
 
-__host__ __device__ cucomp128::cucomp128()
-{
-    real = 0.0;
-    imag = 0.0;
-    return;
-}
-
-__host__ __device__ cucomp128::~cucomp128()
-{
-    real = 0.0;
-    imag = 0.0;
-    return;
-}
-
-__host__ __device__ cucomp128::cucomp128(const cucomp128 &other)
-{
-    real = other.real;
-    imag = other.imag;
-
-    return;
-}
-
-__host__ __device__ cucomp128::cucomp128(const double &other)
-{
-    real = other;
-    imag = 0.0;
-    return;
-}
-
-__host__ __device__ cucomp128& cucomp128::operator=(cucomp128& other)
-{
-    real = other.real;
-    imag = other.imag;
-    return *this;
-}
-
-__host__ __device__ const cucomp128& cucomp128::operator=(const cucomp128& other) 
-{
-    this->real = other.real;
-    this->imag = other.imag;
-    return *this;
-}
-
-__host__ __device__ cucomp128& cucomp128::operator=(double& other)
-{
-    real = other;
-    imag = 0.0;
-    return *this;
-}
-
-__host__ __device__ const cucomp128& cucomp128::operator=(const double& other) 
-{
-    this->real = other;
-    this->imag = 0.0;
-    return *this;
-}
-
-__host__ __device__ double& cucomp128::operator[](int& ind)
-{
-    if(ind==0)
+    __host__ __device__ cucomp128::cucomp128()
     {
-        return this->real;
+        real = 0.0;
+        imag = 0.0;
+        return;
     }
-    else
+    
+    __host__ __device__ cucomp128::~cucomp128()
     {
-        return this->imag;
+        real = 0.0;
+        imag = 0.0;
+        return;
     }
-}
-
-__host__ __device__ const double& cucomp128::operator[](const int& ind) const
-{
-    if(ind==0)
+    
+    __host__ __device__ cucomp128::cucomp128(const cucomp128 &other)
     {
-        return this->real;
+        real = other.real;
+        imag = other.imag;
+    
+        return;
     }
-    else
+    
+    __host__ __device__ cucomp128::cucomp128(const double &other)
     {
-        return this->imag;
+        real = other;
+        imag = 0.0;
+        return;
     }
-}
-
-__host__ __device__ cucomp128 cucomp128::operator+(const cucomp128& z)
-{
-  cucomp128 ret;
-  ret.real = real + z.real;
-  ret.imag = imag + z.imag;
-  return ret;  
-}
-
-__host__ __device__ cucomp128 cucomp128::operator-(const cucomp128& z)
-{
-  cucomp128 ret;
-  ret.real = real - z.real;
-  ret.imag = imag - z.imag;
-  return ret;  
-}
-
-__host__ __device__ cucomp128 cucomp128::operator*(const cucomp128& z)
-{
-  cucomp128 ret;
-  ret.real = (real*z.real - imag*z.imag);
-  ret.imag = (imag*z.real + real*z.imag);
-  return ret;  
-}
-
-__host__ __device__ cucomp128 cucomp128::operator/(const cucomp128& z)
-{
-  cucomp128 ret;
-  double zm2 = z.real*z.real+z.imag*z.imag;
-  
-  if(zm2>0.0)
-  {
-      ret.real = (this->real*z.real+this->imag*z.imag)/zm2;
-      ret.imag = (this->imag*z.real-this->real*z.imag)/zm2;
-  }
-  else
-  {
-      ret.real = (double) finf;
-      ret.imag = (double) finf;
-  }
-
-  return ret;  
-}
-
-__host__ __device__ cucomp128 cucomp128::operator+(const double& z)
-{
-    cucomp128 ret;
-    ret.real = this->real + z;
-    ret.imag = this->imag;
-    return ret;
-}
-__host__ __device__ cucomp128 cucomp128::operator-(const double& z)
-{
-    cucomp128 ret;
-    ret.real = real-z;
-    ret.imag = imag;
-    return ret;
-}
-__host__ __device__ cucomp128 cucomp128::operator*(const double& z)
-{
-    cucomp128 ret;
-    ret.real = real*z;
-    ret.imag = imag*z;
-    return ret;
-}
-__host__ __device__ cucomp128 cucomp128::operator/(const double& z)
-{
-    cucomp128 ret;
-    if(z!=0.0f)
+    
+    __host__ __device__ cucomp128::cucomp128(const double &_real, const double &_imag)
     {
-        ret.real = real/z;
-        ret.imag = imag/z;
+        real = _real;
+        imag = _imag;
+        return;
     }
-    else
+    
+    __host__ __device__ cucomp128& cucomp128::operator=(const cucomp128& other)
     {
-        ret.real = finf;
-        ret.imag = finf;
+        real = other.real;
+        imag = other.imag;
+        return *this;
     }
-
-    return ret;
-}
-
-__host__ __device__ bool cucomp128::operator==(const cucomp128& z) const
-{
-    bool ret = 0;
-    if(z.real == real && z.imag == imag)
+    
+    
+    __host__ __device__ cucomp128& cucomp128::operator=(const double& other)
     {
-        ret = 1;
+        real = other;
+        imag = 0.0;
+        return *this;
     }
-    return ret;
-}
-
-__host__ __device__ bool cucomp128::operator!=(const cucomp128& z) const
-{
-    return !(*this==z);
-}
-
-//sort first by real value, then by imaginary value
-//this is done so that an ordering exists, as long as two values aren't equal
-__host__ __device__ bool cucomp128::operator>(const cucomp128& z) const
-{
-    bool ret = 0;
-    if(this->real>z.real)
+    
+    __host__ __device__ double& cucomp128::operator[](int& ind)
     {
-        ret = 1;
+        if(ind==0)
+        {
+            return this->real;
+        }
+        else
+        {
+            return this->imag;
+        }
     }
-    else if(this->real==z.real)
+    
+    __host__ __device__ const double& cucomp128::operator[](const int& ind) const
     {
-        if(this->imag>z.imag)
+        if(ind==0)
+        {
+            return this->real;
+        }
+        else
+        {
+            return this->imag;
+        }
+    }
+    
+    __host__ __device__ cucomp128 cucomp128::operator+(const cucomp128& z) const
+    {
+      cucomp128 ret;
+      ret.real = real + z.real;
+      ret.imag = imag + z.imag;
+      return ret;  
+    }
+    
+    __host__ __device__ cucomp128 cucomp128::operator-(const cucomp128& z) const
+    {
+      cucomp128 ret;
+      ret.real = real - z.real;
+      ret.imag = imag - z.imag;
+      return ret;  
+    }
+    
+    __host__ __device__ cucomp128 cucomp128::operator*(const cucomp128& z) const
+    {
+      cucomp128 ret;
+      ret.real = real*z.real - imag*z.imag;
+      ret.imag = imag*z.real + real*z.imag;
+      return ret;  
+    }
+    
+    __host__ __device__ cucomp128 cucomp128::operator/(const cucomp128& z) const
+    {
+      cucomp128 ret;
+      double zm2 = z.real*z.real+z.imag*z.imag;
+      
+    //   if(zm2>0.0)
+    //   {
+          ret.real = (real*z.real+imag*z.imag)/zm2;
+          ret.imag = (imag*z.real-real*z.imag)/zm2;
+    //   }
+    //   else
+    //   {
+    //       ret.real = (double) finf;
+    //       ret.imag = (double) finf;
+    //   }
+    
+      return ret;  
+    }
+    
+    __host__ __device__ cucomp128 operator+(const cucomp128& z1, const double& z2)
+    {
+        cucomp128 ret;
+        ret.real = z1.real+z2;
+        ret.imag = z1.imag;
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 operator-(const cucomp128& z1, const double& z2)
+    {
+        cucomp128 ret;
+        ret.real = z1.real-z2;
+        ret.imag = z1.imag;
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 operator*(const cucomp128& z1, const double& z2)
+    {
+        cucomp128 ret;
+        ret.real = z1.real*z2;
+        ret.imag = z1.imag*z2;
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 operator/(const cucomp128& z1, const double& z2)
+    {
+        cucomp128 ret;
+        ret.real = z1.real/z2;
+        ret.imag = z1.imag/z2;
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 operator+(const double& z1, const cucomp128& z2)
+    {
+        cucomp128 ret;
+        ret.real = z1 + z2.real;
+        ret.imag = z2.imag;
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 operator-(const double& z1, const cucomp128& z2)
+    {
+        cucomp128 ret;
+        ret.real = z1 - z2.real;
+        ret.imag = -z2.imag;
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 operator*(const double& z1, const cucomp128& z2)
+    {
+        cucomp128 ret;
+        ret.real = z1*z2.real;
+        ret.imag = z1*z2.imag;
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 operator/(const double& z1, const cucomp128& z2)
+    {
+        cucomp128 ret;
+        double zmg = z2.real*z2.real + z2.imag*z2.imag;
+        ret.real = z1*(z2.real)/zmg;
+        ret.imag = -z1*(z2.imag)/zmg;
+        return ret;
+    }
+    
+    __host__ __device__ bool cucomp128::operator==(const cucomp128& z) const
+    {
+        bool ret = 0;
+        if(z.real == real && z.imag == imag)
         {
             ret = 1;
         }
+        return ret;
+    }
+    
+    __host__ __device__ bool cucomp128::operator!=(const cucomp128& z) const
+    {
+        return !(*this==z);
+    }
+    
+    //sort first by real value, then by imaginary value
+    //this is done so that an ordering exists, as long as two values aren't equal
+    __host__ __device__ bool cucomp128::operator>(const cucomp128& z) const
+    {
+        bool ret = 0;
+        if(this->real>z.real)
+        {
+            ret = 1;
+        }
+        else if(this->real==z.real)
+        {
+            if(this->imag>z.imag)
+            {
+                ret = 1;
+            }
+            else
+            {
+                ret = 0;
+            }
+        }
         else
         {
             ret = 0;
         }
+        return ret;
     }
-    else
+    
+    __host__ __device__ bool cucomp128::operator<(const cucomp128& z) const
     {
-        ret = 0;
-    }
-    return ret;
-}
-
-__host__ __device__ bool cucomp128::operator<(const cucomp128& z) const
-{
-    bool ret = 0;
-    if(this->real<z.real)
-    {
-        ret = 1;
-    }
-    else if(this->real==z.real)
-    {
-        if(this->imag<z.imag)
+        bool ret = 0;
+        if(this->real<z.real)
         {
             ret = 1;
         }
+        else if(this->real==z.real)
+        {
+            if(this->imag<z.imag)
+            {
+                ret = 1;
+            }
+            else
+            {
+                ret = 0;
+            }
+        }
         else
         {
             ret = 0;
         }
+        return ret;
     }
-    else
+    
+    __host__ __device__ cucomp128 operator-(const cucomp128 &z)
     {
-        ret = 0;
+        cucomp128 ret;
+        ret.real = -z.real;
+        ret.imag = -z.imag;
+        return ret;
     }
-    return ret;
-}
-
-__host__ __device__ cucomp128 operator-(const cucomp128 &z)
-{
-	cucomp128 ret;
-	ret.real = -z.real;
-	ret.imag = -z.imag;
-	return ret;
-}
-
-__host__ __device__ bool cucomp128::operator>=(const cucomp128& z) const
-{
-    bool ret = (*this==z || *this>z);
-    return ret;
-}
-
-__host__ __device__ bool cucomp128::operator<=(const cucomp128& z) const
-{
-    bool ret = (*this==z || *this<z);
-    return ret;
-}
-
-__host__ __device__ bool cucomp128::isnan() const
-{
-  bool ret = 0;
-  if(::isnan(this->real) || ::isnan(this->imag))
-  {
-    ret = 1;
-  }
-  return ret;
-}
-
-__host__ __device__ bool cucomp128::isinf() const
-{
-  bool ret = 0;
-  //calls math.h isinf()
-  if(::isinf(this->real) || ::isinf(this->imag))
-  {
-    ret = 1;
-  }
-  return ret;
-}
-
-__host__ __device__ bool cucomp128::isreal() const
-{
-    bool ret = 1;
-    if(imag!=0.0f)
+    
+    __host__ __device__ bool cucomp128::operator>=(const cucomp128& z) const
     {
-        ret = 0;
+        bool ret = (*this==z || *this>z);
+        return ret;
     }
-    return ret;
-}
-
-__host__ __device__ bool cucomp128::isimag() const
-{
-    bool ret = 1;
-    if(real!=0.0f)
+    
+    __host__ __device__ bool cucomp128::operator<=(const cucomp128& z) const
     {
-        ret = 0;
+        bool ret = (*this==z || *this<z);
+        return ret;
     }
-    return ret;
-}
-
-__host__ __device__ bool cucomp128::iszero() const
-{
-    bool ret = 1;
-    if(real!=0.0f || imag!=0.0f)
+    
+    __host__ __device__ bool cucomp128::isnan() const
     {
-        ret = 0;
+      bool ret = 0;
+      if(::isnan(this->real) || ::isnan(this->imag))
+      {
+        ret = 1;
+      }
+      return ret;
+    }
+    
+    __host__ __device__ bool cucomp128::isinf() const
+    {
+      bool ret = 0;
+      //calls math.h isinf()
+      if(::isinf(this->real) || ::isinf(this->imag))
+      {
+        ret = 1;
+      }
+      return ret;
+    }
+    
+    __host__ __device__ bool cucomp128::isreal() const
+    {
+        bool ret = 1;
+        if(imag!=0.0)
+        {
+            ret = 0;
+        }
+        return ret;
+    }
+    
+    __host__ __device__ bool cucomp128::isimag() const
+    {
+        bool ret = 1;
+        if(real!=0.0)
+        {
+            ret = 0;
+        }
+        return ret;
+    }
+    
+    __host__ __device__ bool cucomp128::iszero() const
+    {
+        bool ret = 1;
+        if(real!=0.0 || imag!=0.0)
+        {
+            ret = 0;
+        }
+        return ret;
+    }
+    
+    __host__ __device__ double cucomp128::arg() const
+    {
+        double ret = 0.0;
+        ret = (double) amscuda::arg((double)real,(double)imag);
+        return ret;
+    }
+    
+    __host__ __device__ double cucomp128::mag() const
+    {
+        double ret = 0.0;
+        ret = ::sqrt(real*real+imag*imag);
+        return ret;
+    }
+    
+    __host__ __device__ double cucomp128::abs() const
+    {
+        double ret = 0.0;
+        ret = ::sqrt(real*real+imag*imag);
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 cucomp128::conj() const
+    {
+        cucomp128 ret;
+        ret.real = real;
+        ret.imag = -imag;
+        return ret;
+    }
+    
+    __host__ __device__ double arg(const cucomp128 &z)
+    {
+        return z.arg();
+    }
+    
+    __host__ __device__ double abs(const cucomp128 &z)
+    {
+        return z.mag();
+    }
+    
+    __host__ __device__ double real(const cucomp128 &z)
+    {
+        return z.real;
+    }
+    
+    __host__ __device__ double imag(const cucomp128 &z)
+    {
+        return z.imag;
+    }
+    
+    __host__ __device__ cucomp128 sin(const cucomp128 &z)
+    {
+        cucomp128 ret;
+        cucomp128 im1 = cucomp128(0.0,1.0f);
+        cucomp128 div = cucomp128(0.0,2.0f);
+    
+        ret = (exp(im1*z)-exp(-im1*z))/div;
+    
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 cos(const cucomp128 &z)
+    {
+        cucomp128 ret;
+        cucomp128 im1 = cucomp128(0.0,1.0f);
+        cucomp128 div = cucomp128(2.0f,0.0);
+    
+        ret = (exp(im1*z)+exp(-im1*z))/div;
+    
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 tan(const cucomp128 &z)
+    {
+        return sin(z)/cos(z);
+    }
+    
+    __host__ __device__ cucomp128 exp(const cucomp128 &z)
+    {
+        cucomp128 ret;
+        ret.real = ::exp(z.real)*::cos(z.imag);
+        ret.imag = ::exp(z.real)*::sin(z.imag);
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 log(const cucomp128 &z)
+    {
+        cucomp128 ret;
+        ret.real = ::log(::sqrt(z.real*z.real+z.imag*z.imag));
+        ret.imag = amscuda::arg(z.real,z.imag);
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 conj(const cucomp128 &z)
+    {
+        return z.conj();
+    }
+    
+    __host__ __device__ cucomp128 cosh(const cucomp128 &z)
+    {
+        cucomp128 ret;
+        cucomp128 div = cucomp128(2.0f,0.0);
+    
+        ret = (exp(z)+exp(-z))/div;
+    
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 sinh(const cucomp128 &z)
+    {
+        cucomp128 ret;
+        cucomp128 div = cucomp128(2.0f,0.0);
+    
+        ret = (exp(z)-exp(-z))/div;
+    
+        return ret;
+    }
+    
+    __host__ __device__ cucomp128 tanh(const cucomp128 &z)
+    {
+        return sinh(z)/cosh(z);
+    }
+    
+    __host__ __device__ cucomp128 pow(const cucomp128 &z1, const cucomp128 &z2)
+    {
+        cucomp128 ret;
+        if(z1.mag()>0.0)
+            ret = amscuda::cmp::exp(amscuda::cmp::log(z1)*z2);
+        else
+            ret = cucomp128(0.0,0.0);
+        return ret;
+    }
+    
+    // //returns "complex sign" of complex number - 0, or a unit number with same argument
+    __host__ __device__ cucomp128 csgn(const cucomp128 &z)
+    {
+        cucomp128 ret;
+        double mag = z.mag();
+        if(mag>0.0)
+        {
+            ret = z;
+            ret = ret/mag;
+        }
+        return ret;
+    }
+    
+    //////////////////////////
+    //accumulation operators//
+    //////////////////////////
+    
+    __host__ __device__ cucomp128& cucomp128::operator+=(const cucomp128& z)
+    {
+        real += z.real;
+        imag += z.imag;
+        return (*this);
+    }
+    
+    __host__ __device__ cucomp128& cucomp128::operator-=(const cucomp128& z)
+    {
+        real -= z.real;
+        imag -= z.imag;
+        return (*this);
+    }
+    
+    __host__ __device__ cucomp128& cucomp128::operator*=(const cucomp128& z)
+    {
+        cucomp128 t = (*this);
+        real = t.real*z.real - t.imag*z.imag;
+        imag = t.real*z.imag + t.imag*z.real;
+        return (*this);
+    }
+    
+    __host__ __device__ cucomp128& cucomp128::operator/=(const cucomp128& z)
+    {
+        amscuda::cmp::cucomp128 t = *this;
+        *this = t/z;
+        return (*this);
+    }
+    
+    __host__ __device__ cucomp128& cucomp128::operator+=(const double& z)
+    {
+        real += z;
+        return (*this);
+    }
+    
+    __host__ __device__ cucomp128& cucomp128::operator-=(const double& z)
+    {
+        real -= z;
+        return (*this);
+    }
+    
+    __host__ __device__ cucomp128& cucomp128::operator*=(const double& z)
+    {
+        real *= z;
+        imag *= z;
+        return (*this);
+    }
+    
+    __host__ __device__ cucomp128& cucomp128::operator/=(const double& z)
+    {
+        real /= z;
+        imag /= z;
+        return (*this);
     }
-    return ret;
-}
-
-__host__ __device__ double cucomp128::arg() const
-{
-    double ret = 0.0;
-    ret = (double) amscuda::arg((double)real,(double)imag);
-    return ret;
-}
-
-__host__ __device__ double cucomp128::mag() const
-{
-    double ret = 0.0;
-    ret = ::sqrt(real*real+imag*imag);
-    return ret;
-}
-
-__host__ __device__ cucomp128 cucomp128::conj() const
-{
-    cucomp128 ret;
-    ret.real = real;
-    ret.imag = -imag;
-    return ret;
-}
-
-__host__ __device__ double arg(cucomp128 z)
-{
-    return z.arg();
-}
-
-__host__ __device__ double abs(cucomp128 z)
-{
-    return z.mag();
-}
-
-__host__ __device__ cucomp128 dtocomp(double _r, double _i)
-{
-    cucomp128 ret;
-    ret.real = _r;
-    ret.imag = _i;
-    return ret;
-}
-
-__host__ __device__ double real(cucomp128 z)
-{
-    return z.real;
-}
-
-__host__ __device__ double imag(cucomp128 z)
-{
-    return z.imag;
-}
-
-__host__ __device__ cucomp128 sin(cucomp128 z)
-{
-    cucomp128 ret;
-    cucomp128 im1 = dtocomp(0.0f,1.0f);
-    cucomp128 div = dtocomp(0.0f,2.0f);
-
-    ret = (exp(im1*z)-exp(-im1*z))/div;
-
-    return ret;
-}
-
-__host__ __device__ cucomp128 cos(cucomp128 z)
-{
-    cucomp128 ret;
-    cucomp128 im1 = dtocomp(0.0f,1.0f);
-    cucomp128 div = dtocomp(2.0f,0.0f);
-
-    ret = (exp(im1*z)+exp(-im1*z))/div;
-
-    return ret;
-}
-
-__host__ __device__ cucomp128 tan(cucomp128 z)
-{
-    return sin(z)/cos(z);
-}
-
-__host__ __device__ cucomp128 exp(cucomp128 z)
-{
-    cucomp128 ret;
-    ret.real = ::exp(z.real)*::cos(z.imag);
-    ret.imag = ::exp(z.real)*::sin(z.imag);
-    return ret;
-}
-
-__host__ __device__ cucomp128 log(cucomp128 z)
-{
-    cucomp128 ret;
-    ret.real = ::log(::sqrt(z.real*z.real+z.imag*z.imag));
-    ret.imag = amscuda::arg(z.real,z.imag);
-    return ret;
-}
-
-__host__ __device__ cucomp128 conj(cucomp128 z)
-{
-    return z.conj();
-}
-
-__host__ __device__ cucomp128 cosh(cucomp128 z)
-{
-    cucomp128 ret;
-    cucomp128 div = dtocomp(2.0f,0.0f);
-
-    ret = (exp(z)+exp(-z))/div;
-
-    return ret;
-}
-
-__host__ __device__ cucomp128 sinh(cucomp128 z)
-{
-    cucomp128 ret;
-    cucomp128 div = dtocomp(2.0f,0.0f);
-
-    ret = (exp(z)-exp(-z))/div;
-
-    return ret;
-}
-
-__host__ __device__ cucomp128 tanh(cucomp128 z)
-{
-    return sinh(z)/cosh(z);
-}
-
-__host__ __device__ cucomp128 pow(cucomp128 z1, cucomp128 z2)
-{
-	cucomp128 ret;
-	if(z1.mag()>0.0)
-		ret = exp(z2*log(z1));
-	else
-		ret = dtocomp(0.0f,0.0f);
-	return ret;
-}
 
 
 void test_cucomp128_1()
@@ -444,15 +536,15 @@ void test_cucomp128_1()
     cucomp128 z1;
     cucomp128 a,b,c;
     double d1;
-    double f1;
+    float f1;
 
-    printf("sizeof double=%ld\n",(long)(8*sizeof(f1)));
+    printf("sizeof float=%ld\n",(long)(8*sizeof(f1)));
     printf("sizeof double=%ld\n",(long)(8*sizeof(d1)));
     printf("sizeof complex=%ld\n",(long)(8*sizeof(z1)));
     printf("sizeof cucomp128=%ld\n",(long)(8*sizeof(a)));
 
-    a = dtocomp(1.0,1.0);
-    b = dtocomp(1.0,-1.0);
+    a = cucomp128(1.0,1.0);
+    b = cucomp128(1.0,-1.0);
 
     printf("a=%1.4f + %1.4fi\n",a[0],a[1]);
     printf("b=%1.4f + %1.4fi\n",b[0],b[1]);
diff --git a/src/amsculib3/math/amscu_comp64.cu b/src/amsculib3/math/amscu_comp64.cu
index ff848f9..3173eac 100644
--- a/src/amsculib3/math/amscu_comp64.cu
+++ b/src/amsculib3/math/amscu_comp64.cu
@@ -329,7 +329,7 @@ __host__ __device__ bool cucomp64::iszero() const
 __host__ __device__ float cucomp64::arg() const
 {
     float ret = 0.0f;
-    ret = (float) amscuda::arg((double)real,(double)imag);
+    ret = (float) amscuda::argf(real,imag);
     return ret;
 }