diff --git a/numexpr/complex_functions.hpp b/numexpr/complex_functions.hpp
index 6d48866..ae89167 100644
--- a/numexpr/complex_functions.hpp
+++ b/numexpr/complex_functions.hpp
@@ -10,16 +10,17 @@
   See LICENSE.txt for details about copyright and rights to use.
 **********************************************************************/
 
-// TODO: Could just use std::complex<float> and std::complex<double>
+// Replace npy_cdouble with std::complex<double>
+#include <complex>
 
 /* constants */
-static npy_cdouble nc_1 = {1., 0.};
-static npy_cdouble nc_half = {0.5, 0.};
-static npy_cdouble nc_i = {0., 1.};
-static npy_cdouble nc_i2 = {0., 0.5};
+static std::complex<double> nc_1(1., 0.);
+static std::complex<double> nc_half(0.5, 0.);
+static std::complex<double> nc_i(0., 1.);
+static std::complex<double> nc_i2(0., 0.5);
 /*
-static npy_cdouble nc_mi = {0., -1.};
-static npy_cdouble nc_pi2 = {M_PI/2., 0.};
+static std::complex<double> nc_mi = {0., -1.};
+static std::complex<double> nc_pi2 = {M_PI/2., 0.};
 */
 
 /* *************************** WARNING *****************************
@@ -31,42 +32,42 @@ after the first *r has been overwritten.
 */
 
 static void
-nc_assign(npy_cdouble *x, npy_cdouble *r)
+nc_assign(std::complex<double> *x, std::complex<double> *r)
 {
-  r->real = x->real;
-  r->imag = x->imag;
+  r->real(x->real());
+  r->imag(x->imag());
   return;
 }
 
 static void
-nc_sum(npy_cdouble *a, npy_cdouble *b, npy_cdouble *r)
+nc_sum(std::complex<double> *a, std::complex<double> *b, std::complex<double> *r)
 {
-    r->real = a->real + b->real;
-    r->imag = a->imag + b->imag;
+    r->real(a->real() + b->real());
+    r->imag(a->imag() + b->imag());
     return;
 }
 
 static void
-nc_diff(npy_cdouble *a, npy_cdouble *b, npy_cdouble *r)
+nc_diff(std::complex<double> *a, std::complex<double> *b, std::complex<double> *r)
 {
-    r->real = a->real - b->real;
-    r->imag = a->imag - b->imag;
+    r->real(a->real() - b->real());
+    r->imag(a->imag() - b->imag());
     return;
 }
 
 static void
-nc_neg(npy_cdouble *a, npy_cdouble *r)
+nc_neg(std::complex<double> *a, std::complex<double> *r)
 {
-    r->real = -a->real;
-    r->imag = -a->imag;
+    r->real(-a->real());
+    r->imag(-a->imag());
     return;
 }
 
 static void
-nc_conj(npy_cdouble *a, npy_cdouble *r)
+nc_conj(std::complex<double> *a, std::complex<double> *r)
 {
-    r->real = a->real;
-    r->imag = -a->imag;
+    r->real(a->real());
+    r->imag(-a->imag());
     return;
 }
 
@@ -85,109 +86,109 @@ inline double fconj(double x)
 }
 
 static void
-nc_prod(npy_cdouble *a, npy_cdouble *b, npy_cdouble *r)
+nc_prod(std::complex<double> *a, std::complex<double> *b, std::complex<double> *r)
 {
-    double ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
-    r->real = ar*br - ai*bi;
-    r->imag = ar*bi + ai*br;
+    double ar=a->real(), br=b->real(), ai=a->imag(), bi=b->imag();
+    r->real(ar*br - ai*bi);
+    r->imag(ar*bi + ai*br);
     return;
 }
 
 static void
-nc_quot(npy_cdouble *a, npy_cdouble *b, npy_cdouble *r)
+nc_quot(std::complex<double> *a, std::complex<double> *b, std::complex<double> *r)
 {
-    double ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
+    double ar=a->real(), br=b->real(), ai=a->imag(), bi=b->imag();
     double d = br*br + bi*bi;
-    r->real = (ar*br + ai*bi)/d;
-    r->imag = (ai*br - ar*bi)/d;
+    r->real((ar*br + ai*bi)/d);
+    r->imag((ai*br - ar*bi)/d);
     return;
 }
 
 static void
-nc_sqrt(npy_cdouble *x, npy_cdouble *r)
+nc_sqrt(std::complex<double> *x, std::complex<double> *r)
 {
     double s,d;
-    if (x->real == 0. && x->imag == 0.)
+    if (x->real() == 0. && x->imag() == 0.)
         *r = *x;
     else {
-        s = sqrt((fabs(x->real) + hypot(x->real,x->imag))/2);
-        d = x->imag/(2*s);
-        if (x->real > 0.) {
-            r->real = s;
-            r->imag = d;
+        s = sqrt((fabs(x->real()) + hypot(x->real(),x->imag()))/2);
+        d = x->imag()/(2*s);
+        if (x->real() > 0.) {
+            r->real(s);
+            r->imag(d);
         }
-        else if (x->imag >= 0.) {
-            r->real = d;
-            r->imag = s;
+        else if (x->imag() >= 0.) {
+            r->real(d);
+            r->imag(s);
         }
         else {
-            r->real = -d;
-            r->imag = -s;
+            r->real(-d);
+            r->imag(-s);
         }
     }
     return;
 }
 
 static void
-nc_log(npy_cdouble *x, npy_cdouble *r)
+nc_log(std::complex<double> *x, std::complex<double> *r)
 {
-    double l = hypot(x->real,x->imag);
-    r->imag = atan2(x->imag, x->real);
-    r->real = log(l);
+    double l = hypot(x->real(),x->imag());
+    r->imag(atan2(x->imag(), x->real()));
+    r->real(log(l));
     return;
 }
 
 static void
-nc_log1p(npy_cdouble *x, npy_cdouble *r)
+nc_log1p(std::complex<double> *x, std::complex<double> *r)
 {
-    double l = hypot(x->real + 1.0,x->imag);
-    r->imag = atan2(x->imag, x->real + 1.0);
-    r->real = log(l);
+    double l = hypot(x->real() + 1.0,x->imag());
+    r->imag(atan2(x->imag(), x->real() + 1.0));
+    r->real(log(l));
     return;
 }
 
 static void
-nc_exp(npy_cdouble *x, npy_cdouble *r)
+nc_exp(std::complex<double> *x, std::complex<double> *r)
 {
-    double a = exp(x->real);
-    r->real = a*cos(x->imag);
-    r->imag = a*sin(x->imag);
+    double a = exp(x->real());
+    r->real(a*cos(x->imag()));
+    r->imag(a*sin(x->imag()));
     return;
 }
 
 static void
-nc_expm1(npy_cdouble *x, npy_cdouble *r)
+nc_expm1(std::complex<double> *x, std::complex<double> *r)
 {
-    double a = sin(x->imag / 2);
-    double b = exp(x->real);
-    r->real = expm1(x->real) * cos(x->imag) - 2 * a * a;
-    r->imag = b * sin(x->imag);
+    double a = sin(x->imag() / 2);
+    double b = exp(x->real());
+    r->real(expm1(x->real()) * cos(x->imag()) - 2 * a * a);
+    r->imag(b * sin(x->imag()));
     return;
 }
 
 static void
-nc_pow(npy_cdouble *a, npy_cdouble *b, npy_cdouble *r)
+nc_pow(std::complex<double> *a, std::complex<double> *b, std::complex<double> *r)
 {
     npy_intp n;
-    double ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
+    double ar=a->real(), br=b->real(), ai=a->imag(), bi=b->imag();
 
     if (br == 0. && bi == 0.) {
-        r->real = 1.;
-        r->imag = 0.;
+        r->real(1.);
+        r->imag(0.);
         return;
     }
     if (ar == 0. && ai == 0.) {
-        r->real = 0.;
-        r->imag = 0.;
+        r->real(0.);
+        r->imag(0.);
         return;
     }
     if (bi == 0 && (n=(npy_intp)br) == br) {
         if (n > -100 && n < 100) {
-        npy_cdouble p, aa;
+        std::complex<double> p, aa;
         npy_intp mask = 1;
         if (n < 0) n = -n;
         aa = nc_1;
-        p.real = ar; p.imag = ai;
+        p.real(ar); p.imag(ai);
         while (1) {
             if (n & mask)
                 nc_prod(&aa,&p,&aa);
@@ -195,7 +196,7 @@ nc_pow(npy_cdouble *a, npy_cdouble *b, npy_cdouble *r)
             if (n < mask || mask <= 0) break;
             nc_prod(&p,&p,&p);
         }
-        r->real = aa.real; r->imag = aa.imag;
+        r->real(aa.real()); r->imag(aa.imag());
         if (br < 0) nc_quot(&nc_1, r, r);
         return;
         }
@@ -210,19 +211,19 @@ nc_pow(npy_cdouble *a, npy_cdouble *b, npy_cdouble *r)
 
 
 static void
-nc_prodi(npy_cdouble *x, npy_cdouble *r)
+nc_prodi(std::complex<double> *x, std::complex<double> *r)
 {
-    double xr = x->real;
-    r->real = -x->imag;
-    r->imag = xr;
+    double xr = x->real();
+    r->real(-x->imag());
+    r->imag(xr);
     return;
 }
 
 
 static void
-nc_acos(npy_cdouble *x, npy_cdouble *r)
+nc_acos(std::complex<double> *x, std::complex<double> *r)
 {
-    npy_cdouble a, *pa=&a;
+    std::complex<double> a, *pa=&a;
 
     nc_assign(x, pa);
     nc_prod(x,x,r);
@@ -240,9 +241,9 @@ nc_acos(npy_cdouble *x, npy_cdouble *r)
 }
 
 static void
-nc_acosh(npy_cdouble *x, npy_cdouble *r)
+nc_acosh(std::complex<double> *x, std::complex<double> *r)
 {
-    npy_cdouble t, a, *pa=&a;
+    std::complex<double> t, a, *pa=&a;
 
     nc_assign(x, pa);
     nc_sum(x, &nc_1, &t);
@@ -260,9 +261,9 @@ nc_acosh(npy_cdouble *x, npy_cdouble *r)
 }
 
 static void
-nc_asin(npy_cdouble *x, npy_cdouble *r)
+nc_asin(std::complex<double> *x, std::complex<double> *r)
 {
-    npy_cdouble a, *pa=&a;
+    std::complex<double> a, *pa=&a;
     nc_prodi(x, pa);
     nc_prod(x, x, r);
     nc_diff(&nc_1, r, r);
@@ -280,9 +281,9 @@ nc_asin(npy_cdouble *x, npy_cdouble *r)
 
 
 static void
-nc_asinh(npy_cdouble *x, npy_cdouble *r)
+nc_asinh(std::complex<double> *x, std::complex<double> *r)
 {
-    npy_cdouble a, *pa=&a;
+    std::complex<double> a, *pa=&a;
     nc_assign(x, pa);
     nc_prod(x, x, r);
     nc_sum(&nc_1, r, r);
@@ -296,9 +297,9 @@ nc_asinh(npy_cdouble *x, npy_cdouble *r)
 }
 
 static void
-nc_atan(npy_cdouble *x, npy_cdouble *r)
+nc_atan(std::complex<double> *x, std::complex<double> *r)
 {
-    npy_cdouble a, *pa=&a;
+    std::complex<double> a, *pa=&a;
     nc_diff(&nc_i, x, pa);
     nc_sum(&nc_i, x, r);
     nc_quot(r, pa, r);
@@ -311,9 +312,9 @@ nc_atan(npy_cdouble *x, npy_cdouble *r)
 }
 
 static void
-nc_atanh(npy_cdouble *x, npy_cdouble *r)
+nc_atanh(std::complex<double> *x, std::complex<double> *r)
 {
-    npy_cdouble a, b, *pa=&a, *pb=&b;
+    std::complex<double> a, b, *pa=&a, *pb=&b;
     nc_assign(x, pa);
     nc_diff(&nc_1, pa, r);
     nc_sum(&nc_1, pa, pb);
@@ -327,20 +328,20 @@ nc_atanh(npy_cdouble *x, npy_cdouble *r)
 }
 
 static void
-nc_cos(npy_cdouble *x, npy_cdouble *r)
+nc_cos(std::complex<double> *x, std::complex<double> *r)
 {
-    double xr=x->real, xi=x->imag;
-    r->real = cos(xr)*cosh(xi);
-    r->imag = -sin(xr)*sinh(xi);
+    double xr=x->real(), xi=x->imag();
+    r->real(cos(xr)*cosh(xi));
+    r->imag(-sin(xr)*sinh(xi));
     return;
 }
 
 static void
-nc_cosh(npy_cdouble *x, npy_cdouble *r)
+nc_cosh(std::complex<double> *x, std::complex<double> *r)
 {
-    double xr=x->real, xi=x->imag;
-    r->real = cos(xi)*cosh(xr);
-    r->imag = sin(xi)*sinh(xr);
+    double xr=x->real(), xi=x->imag();
+    r->real(cos(xi)*cosh(xr));
+    r->imag(sin(xi)*sinh(xr));
     return;
 }
 
@@ -348,39 +349,39 @@ nc_cosh(npy_cdouble *x, npy_cdouble *r)
 #define M_LOG10_E 0.434294481903251827651128918916605082294397
 
 static void
-nc_log10(npy_cdouble *x, npy_cdouble *r)
+nc_log10(std::complex<double> *x, std::complex<double> *r)
 {
     nc_log(x, r);
-    r->real *= M_LOG10_E;
-    r->imag *= M_LOG10_E;
+    r->real(r->real() * M_LOG10_E);
+    r->imag(r->imag() * M_LOG10_E);
     return;
 }
 
 static void
-nc_sin(npy_cdouble *x, npy_cdouble *r)
+nc_sin(std::complex<double> *x, std::complex<double> *r)
 {
-    double xr=x->real, xi=x->imag;
-    r->real = sin(xr)*cosh(xi);
-    r->imag = cos(xr)*sinh(xi);
+    double xr=x->real(), xi=x->imag();
+    r->real(sin(xr)*cosh(xi));
+    r->imag(cos(xr)*sinh(xi));
     return;
 }
 
 static void
-nc_sinh(npy_cdouble *x, npy_cdouble *r)
+nc_sinh(std::complex<double> *x, std::complex<double> *r)
 {
-    double xr=x->real, xi=x->imag;
-    r->real = cos(xi)*sinh(xr);
-    r->imag = sin(xi)*cosh(xr);
+    double xr=x->real(), xi=x->imag();
+    r->real(cos(xi)*sinh(xr));
+    r->imag(sin(xi)*cosh(xr));
     return;
 }
 
 static void
-nc_tan(npy_cdouble *x, npy_cdouble *r)
+nc_tan(std::complex<double> *x, std::complex<double> *r)
 {
     double sr,cr,shi,chi;
     double rs,is,rc,ic;
     double d;
-    double xr=x->real, xi=x->imag;
+    double xr=x->real(), xi=x->imag();
     sr = sin(xr);
     cr = cos(xr);
     shi = sinh(xi);
@@ -390,18 +391,18 @@ nc_tan(npy_cdouble *x, npy_cdouble *r)
     rc = cr*chi;
     ic = -sr*shi;
     d = rc*rc + ic*ic;
-    r->real = (rs*rc+is*ic)/d;
-    r->imag = (is*rc-rs*ic)/d;
+    r->real((rs*rc+is*ic)/d);
+    r->imag((is*rc-rs*ic)/d);
     return;
 }
 
 static void
-nc_tanh(npy_cdouble *x, npy_cdouble *r)
+nc_tanh(std::complex<double> *x, std::complex<double> *r)
 {
     double si,ci,shr,chr;
     double rs,is,rc,ic;
     double d;
-    double xr=x->real, xi=x->imag;
+    double xr=x->real(), xi=x->imag();
     si = sin(xi);
     ci = cos(xi);
     shr = sinh(xr);
@@ -411,16 +412,16 @@ nc_tanh(npy_cdouble *x, npy_cdouble *r)
     rc = ci*chr;
     ic = si*shr;
     d = rc*rc + ic*ic;
-    r->real = (rs*rc+is*ic)/d;
-    r->imag = (is*rc-rs*ic)/d;
+    r->real((rs*rc+is*ic)/d);
+    r->imag((is*rc-rs*ic)/d);
     return;
 }
 
 static void
-nc_abs(npy_cdouble *x, npy_cdouble *r)
+nc_abs(std::complex<double> *x, std::complex<double> *r)
 {
-    r->real = sqrt(x->real*x->real + x->imag*x->imag);
-    r->imag = 0;
+    r->real(sqrt(x->real()*x->real() + x->imag()*x->imag()));
+    r->imag(0);
 }
 
 #endif // NUMEXPR_COMPLEX_FUNCTIONS_HPP
diff --git a/numexpr/interp_body.cpp b/numexpr/interp_body.cpp
index 87bed7c..09b9da9 100644
--- a/numexpr/interp_body.cpp
+++ b/numexpr/interp_body.cpp
@@ -199,7 +199,7 @@
         #define s3    ((char   *)x3+j*sb3)
         /* Some temporaries */
         double da, db;
-        npy_cdouble ca, cb;
+        std::complex<double> ca, cb;
 
         switch (op) {
 
@@ -432,19 +432,19 @@
                                                 (const MKL_Complex16*)x1,
                                                 (MKL_Complex16*)dest));
 #else
-            VEC_ARG1(ca.real = c1r;
-                     ca.imag = c1i;
+            VEC_ARG1(ca.real(c1r);
+                     ca.imag(c1i);
                      functions_cc[arg2](&ca, &ca);
-                     cr_dest = ca.real;
-                     ci_dest = ca.imag);
+                     cr_dest = ca.real();
+                     ci_dest = ca.imag());
 #endif
-        case OP_FUNC_CCCN: VEC_ARG2(ca.real = c1r;
-                                    ca.imag = c1i;
-                                    cb.real = c2r;
-                                    cb.imag = c2i;
+        case OP_FUNC_CCCN: VEC_ARG2(ca.real(c1r);
+                                    ca.imag(c1i);
+                                    cb.real(c2r);
+                                    cb.imag(c2i);
                                     functions_ccc[arg3](&ca, &cb, &ca);
-                                    cr_dest = ca.real;
-                                    ci_dest = ca.imag);
+                                    cr_dest = ca.real();
+                                    ci_dest = ca.imag());
 
         case OP_REAL_DC: VEC_ARG1(d_dest = c1r);
         case OP_IMAG_DC: VEC_ARG1(d_dest = c1i);
diff --git a/numexpr/interpreter.cpp b/numexpr/interpreter.cpp
index f5c0c94..b106a6a 100644
--- a/numexpr/interpreter.cpp
+++ b/numexpr/interpreter.cpp
@@ -246,7 +246,7 @@ FuncDDDPtr_vml functions_ddd_vml[] = {
 
 
 
-typedef void (*FuncCCPtr)(npy_cdouble*, npy_cdouble*);
+typedef void (*FuncCCPtr)(std::complex<double>*, std::complex<double>*);
 
 FuncCCPtr functions_cc[] = {
 #define FUNC_CC(fop, s, f, ...) f,
@@ -295,7 +295,7 @@ FuncCCPtr_vml functions_cc_vml[] = {
 #endif
 
 
-typedef void (*FuncCCCPtr)(npy_cdouble*, npy_cdouble*, npy_cdouble*);
+typedef void (*FuncCCCPtr)(std::complex<double>*, std::complex<double>*, std::complex<double>*);
 
 FuncCCCPtr functions_ccc[] = {
 #define FUNC_CCC(fop, s, f) f,