 #include <complex>
 
 #include <cusp/verify.h>
+
+#define ASSERT_COMPLEX_ALMOST_EQUAL(X,Y) {unittest::assert_almost_equal((X.real()),(Y.real()), __FILE__, __LINE__);unittest::assert_almost_equal((X.imag()),(Y.imag()), __FILE__, __LINE__);}
+
+#ifdef __GNUC__
+template< typename T1, typename T2 >
+struct is_same_type      { enum { result = false }; };
+
+template< typename T>
+struct is_same_type<T,T> { enum { result = true }; };
+
+extern "C"{
+  float __complex__ cacosf(float __complex__ z);
+  double __complex__ cacos(double __complex__ z);
+  float __complex__ casinf(float __complex__ z);
+  double __complex__ casin(double __complex__ z);
+  float __complex__ catanf(float __complex__ z);
+  double __complex__ catan(double __complex__ z);
+  float __complex__ cacoshf(float __complex__ z);
+  double __complex__ cacosh(double __complex__ z);
+  float __complex__ casinhf(float __complex__ z);
+  double __complex__ casinh(double __complex__ z);
+  float __complex__ catanhf(float __complex__ z);
+  double __complex__ catanh(double __complex__ z);
+  double creal(double __complex__ z);
+  double cimag(double __complex__ z);
+}
+#endif
+
 
 // Macro to create host and device versions of a unit test
 #define DECLARE_NUMERIC_UNITTEST(VTEST)                    \

  cudaMemcpy(&a,d_a,sizeof(cusp::complex<ValueType>),cudaMemcpyDeviceToHost);
  std::complex<ValueType> b(a.real(),a.imag());
  a = test_complex_compilation_entry<ValueType>();
  ASSERT_ALMOST_EQUAL(a.real(),b.real());
}
DECLARE_NUMERIC_UNITTEST(test_complex_compilation);



void VTEST##Double(void)   {  VTEST<double>(); }  \
DECLARE_UNITTEST(VTEST##Float);                                 \
DECLARE_UNITTEST(VTEST##Double);

template <typename ValueType>
__host__ bool compareWithStd(cusp::complex<ValueType> a){
  //  std::cout << "Testing " << a << std::endl; 
  cusp::complex<ValueType> b(a.real(),a.imag());
  std::complex<ValueType> s_a(a.real(),a.imag());
  std::complex<ValueType> s_b(b.real(),b.imag());
  ASSERT_COMPLEX_ALMOST_EQUAL(a,s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b += a;
  s_b += s_a;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b -= a;
  s_b -= s_a;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b *= a;
  s_b *= s_a;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b /= a;
  s_b /= s_a;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

  b = a*ValueType(2.0);
  s_b = s_a*ValueType(2.0);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = ValueType(2.0)*a;
  s_b = ValueType(2.0)*s_a;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = b*a;
  s_b = s_b*s_a;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = a/ValueType(2.0);
  s_b = s_a/ValueType(2.0);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = ValueType(2.0)/a;
  s_b = ValueType(2.0)/s_a;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = a/b;
  s_b = s_a/s_b;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = a-ValueType(2.0);
  s_b = s_a-ValueType(2.0);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = ValueType(2.0)-a;
  s_b = ValueType(2.0)-s_a;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = a-b;
  s_b = s_a-s_b;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = a+ValueType(2.0);
  s_b = s_a+ValueType(2.0);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = ValueType(2.0)+a;
  s_b = ValueType(2.0)+s_a;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = a+b;
  s_b = s_a+s_b;
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

  b = cusp::abs(a);
  s_b = std::abs(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::arg(a);
  s_b = std::arg(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::norm(a);
  s_b = std::norm(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::conj(a);
  s_b = std::conj(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::polar(norm(a),a.imag());
  s_b = std::polar(norm(a),a.imag());
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

  b = cusp::cos(a);
  s_b = std::cos(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::cosh(a);
  s_b = std::cosh(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::exp(a);
  s_b = std::exp(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

  b = cusp::log(a);
  s_b = std::log(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::log10(a);
  s_b = std::log10(s_a);
  // There's a difference here but the cusp value came closer
  // to the value calculated with arbitrary precision
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

  b = cusp::pow(a,b);
  s_b = std::pow(s_a,s_b);
  // There's a difference here but the cusp value came closer
  // to the value calculated with arbitrary precision
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::pow(a,ValueType(1.3));
  s_b = std::pow(s_a,ValueType(1.3));
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::pow(ValueType(1.4),a);
  s_b = std::pow(ValueType(1.4),s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  // Missing int implentation in std
  b = cusp::pow(a,4);
  s_b = std::pow(s_a,ValueType(4.0));
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::sin(a);
  s_b = std::sin(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::sinh(a);
  s_b = std::sinh(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::sqrt(a);
  s_b = std::sqrt(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::tan(a);  
  s_b = std::tan(s_a);
  // There's a difference here and it's a close call 
  // but I think the cusp value came closer
  // to the value calculated with arbitrary precision
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::tanh(a);
  s_b = std::tanh(s_a);
  // Here we're both clearly the wrong
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

  // Inverse trigonometic functions not part of standard 
  /*
  a = cusp::acos(b);
  s_a = std::acos(s_b);
  ASSERT_COMPLEX_ALMOST_EQUAL(a,s_a);
  b = cusp::asin(a);
  s_b = std::asin(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  a = cusp::atan(b);
  s_a = std::atan(s_b);
  ASSERT_COMPLEX_ALMOST_EQUAL(a,s_a);
  */

#ifdef __GNUC__  
  /* Use the c99 complex function from gcc to test the
   function not part of the standard */
  if(sizeof(ValueType) == 4){
    __complex__ float g_a;
    __complex__ float g_b;
    g_a = s_a.real() + s_a.imag()*__I__;
    g_b = s_b.real() + s_b.imag()*__I__;
    // Ok here we screwed big time
    b = cusp::acos(a);
    g_b = cacosf(g_a);
    s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
    ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

    g_a = s_a.real() + s_a.imag()*__I__;
    g_b = s_b.real() + s_b.imag()*__I__;
    // Ok here we screwed big time
    b = cusp::asin(a);
    g_b = casinf(g_a);
    s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
    ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

    g_a = s_a.real() + s_a.imag()*__I__;
    g_b = s_b.real() + s_b.imag()*__I__;
    // Ok here we screwed big time
    b = cusp::atan(a);
    g_b = catanf(g_a);
    s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
    ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

    g_a = s_a.real() + s_a.imag()*__I__;
    g_b = s_b.real() + s_b.imag()*__I__;
    // Ok here we screwed big time
    b = cusp::acosh(a);
    g_b = cacoshf(g_a);
    s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
    ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

    g_a = s_a.real() + s_a.imag()*__I__;
    g_b = s_b.real() + s_b.imag()*__I__;
    b = cusp::asinh(a);
    g_b = casinhf(g_a);
    s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
    ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

    if(a != cusp::complex<ValueType>(1,0)){
      g_a = s_a.real() + s_a.imag()*__I__;
      g_b = s_b.real() + s_b.imag()*__I__;
      b = cusp::atanh(a);
      g_b = catanhf(g_a);
      s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
      ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
    }
  }else if(sizeof(ValueType) == 8){
    __complex__ double g_a;
    __complex__ double g_b;
    g_a = s_a.real() + s_a.imag()*__I__;
    g_b = s_b.real() + s_b.imag()*__I__;
    // Ok here we screwed big time
    b = cusp::acos(a);
    g_b = cacos(g_a);
    s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
    ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

    g_a = s_a.real() + s_a.imag()*__I__;
    g_b = s_b.real() + s_b.imag()*__I__;
    // Ok here we screwed big time
    b = cusp::asin(a);
    g_b = casin(g_a);
    s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
    ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

    g_a = s_a.real() + s_a.imag()*__I__;
    g_b = s_b.real() + s_b.imag()*__I__;
    // Ok here we screwed big time
    b = cusp::atan(a);
    g_b = catan(g_a);
    s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
    ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

    g_a = s_a.real() + s_a.imag()*__I__;
    g_b = s_b.real() + s_b.imag()*__I__;
    // Ok here we screwed big time
    b = cusp::acosh(a);
    g_b = cacosh(g_a);
    s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
    ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

    g_a = s_a.real() + s_a.imag()*__I__;
    g_b = s_b.real() + s_b.imag()*__I__;
    b = cusp::asinh(a);
    g_b = casinh(g_a);
    s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
    ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

    if(a != cusp::complex<ValueType>(1,0)){
      g_a = s_a.real() + s_a.imag()*__I__;
      g_b = s_b.real() + s_b.imag()*__I__;
      b = cusp::atanh(a);
      g_b = catanh(g_a);
      s_b = std::complex<ValueType>(creal(g_b),cimag(g_b));
      ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
    }
  }
#endif
  return true;
}

template <typename ValueType>
__host__ __device__ cusp::complex<ValueType> test_complex_members(){
