 
     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));

    __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);
  b = cusp::tan(a);  
  s_b = std::tan(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);
  b = cusp::tanh(a);
  s_b = std::tanh(s_a);
  ASSERT_COMPLEX_ALMOST_EQUAL(b,s_b);

  // Inverse trigonometic functions not part of standard 
