diff -uN libcerf-1.3/lib/defs.h.ORG libcerf-1.3/lib/defs.h
--- libcerf-1.3/lib/defs.h.ORG	Fri Feb  1 18:00:14 2013
+++ libcerf-1.3/lib/defs.h	Fri Jul 19 16:49:27 2013
@@ -20,9 +20,13 @@
  * Website:
  *   http://apps.jcns.fz-juelich.de/libcerf
  */
+#ifndef __LIBCERF_DEFS_H__
+#define __LIBCERF_DEFS_H__
 
 
+#ifndef __sun__
 typedef double _Complex cmplx;
+#endif
 
 /* Constructing complex numbers like 0+i*NaN is problematic in C99
    without the C11 CMPLX macro, because 0.+I*NAN may give NaN+i*NAN if
@@ -36,11 +40,21 @@
    glibc and compile with -std=c11... note that icc lies about being
    gcc and probably doesn't have this builtin(?), so exclude icc explicitly */
 
+#ifndef __sun
 #if !defined(CMPLX) && \
     ( __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) && \
     !(defined(__ICC) || defined(__INTEL_COMPILER))
 #  define CMPLX(a,b) __builtin_complex((double) (a), (double) (b))
 #endif
+#else /* __sun__ */
+# define CMPLX(a,b) (cmplx){ {(a),(b)} }
+# if 0
+#  include <values.h>
+#  define INFINITY MAXDOUBLE
+# else
+#  define INFINITY (1./0.)
+# endif
+#endif /* __sun__ */
 
 #ifdef CMPLX // C11
 #  define C(a,b) CMPLX(a,b)
@@ -55,3 +69,5 @@
 #  define Inf (1./0.) 
 #  define NaN (0./0.) 
 #endif
+
+#endif /* __LIBCERF_DEFS_H__ */
diff -uN libcerf-1.3/lib/cerf.h.ORG libcerf-1.3/lib/cerf.h
--- libcerf-1.3/lib/cerf.h.ORG	Mon Feb  4 22:04:18 2013
+++ libcerf-1.3/lib/cerf.h	Fri Jul 19 16:49:43 2013
@@ -43,7 +43,9 @@
  *   w_of_z(3), dawson(3), voigt(3), cerf(3), erfcx(3), erfi(3)
  */
 
+#ifndef __sun__
 #include <complex.h> // C99 complex-number support
+#endif
 
 #ifndef __CERF_H
 #  define __CERF_H
@@ -58,28 +60,47 @@
 #  endif
 __BEGIN_DECLS
 
+#ifndef __sun__
+# define DCMPLX   double _Complex
+#else
+typedef struct { double dat[2]; } cmplx;
+# define DCMPLX   cmplx
+#endif
+
+#ifdef __sun__
+# define creal(z) ((z).dat[0])
+# define cimag(z) ((z).dat[1])
+# include <ieeefp.h>
+#if 1
+# define isinf(x) (!finite(x) && (x)==(x))
+#else
+# define isinf(x) (!finite(x) && !isnand(x))
+# /* no change for test */
+#endif
+#endif
+
 // compute w(z) = exp(-z^2) erfc(-iz), Faddeeva's scaled complex error function
-double _Complex w_of_z   (double _Complex z);
+DCMPLX w_of_z   (DCMPLX z);
 double         im_w_of_x(double x); // special case Im[w(x)] of real x
 double re_w_of_z( double x, double y );
 double im_w_of_z( double x, double y );
 
 // compute erf(z), the error function of complex arguments
-double _Complex cerf  (double _Complex z);
+DCMPLX cerf  (DCMPLX z);
 
 // compute erfc(z) = 1 - erf(z), the complementary error function
-double _Complex cerfc(double _Complex z);
+DCMPLX cerfc(DCMPLX z);
 
 // compute erfcx(z) = exp(z^2) erfc(z), an underflow-compensated version of erfc
-double _Complex cerfcx(double _Complex z);
+DCMPLX cerfcx(DCMPLX z);
 double         erfcx (double x); // special case for real x
 
 // compute erfi(z) = -i erf(iz), the imaginary error function
-double _Complex cerfi(double _Complex z);
+DCMPLX cerfi(DCMPLX z);
 double         erfi (double x); // special case for real x
 
 // compute dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z), Dawson's integral
-double _Complex cdawson(double _Complex z);
+DCMPLX cdawson(DCMPLX z);
 double         dawson (double x); // special case for real x
 
 // compute voigt(x,??), the convolution of a Gaussian and a Lorentzian
diff -uN libcerf-1.3/lib/err_fcts.c.ORG libcerf-1.3/lib/err_fcts.c
--- libcerf-1.3/lib/err_fcts.c.ORG	Mon Feb  4 22:04:18 2013
+++ libcerf-1.3/lib/err_fcts.c	Fri Jul 19 16:50:45 2013
@@ -152,7 +152,11 @@
             return gam / pi / (x*x + gam*gam);
         } else {
             // Regular case, both parameters are nonzero
+#ifndef __sun__
             cmplx z = C(x,gam) / sqrt(2) / sig;
+#else
+            cmplx z = C(x/sqrt(2)/sig, gam/sqrt(2)/sig);
+#endif
             return creal( w_of_z(z) ) / s2pi / sig;
         }
     }
@@ -185,7 +189,11 @@
     double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
     double mIm_z2 = -2*x*y; // Im(-z^2)
     if (mRe_z2 < -750) // underflow
+#ifndef __sun__
         return (x >= 0 ? 1.0 : -1.0);
+#else
+        return (x >= 0 ? C(1.0, 0.0) : C(-1.0, 0.0));
+#endif
 
     /* Handle positive and negative x via different formulas,
        using the mirror symmetries of w, to avoid overflow/underflow
@@ -199,9 +207,18 @@
         }
         /* don't use complex exp function, since that will produce spurious NaN
            values when multiplying w in an overflow situation. */
+#ifndef __sun__
         return 1.0 - exp(mRe_z2) *
             (C(cos(mIm_z2), sin(mIm_z2))
              * w_of_z(C(-y,x)));
+#else
+	{
+	    cmplx zz = w_of_z(C(-y,x));
+	    double RR = exp(mRe_z2), Cv = cos(mIm_z2), Sv = sin(mIm_z2);
+	    return C(1.0 - RR*(Cv*creal(zz) - Sv*cimag(zz)),
+		     - RR*(Cv*cimag(zz) + Sv*creal(zz)));
+	}
+#endif
     }
     else { // x < 0
         if (x > -8e-2) { // duplicate from above to avoid fabs(x) call
@@ -214,21 +231,54 @@
             return C(NaN, y == 0 ? 0 : NaN);
         /* don't use complex exp function, since that will produce spurious NaN
            values when multiplying w in an overflow situation. */
+#ifndef __sun__
         return exp(mRe_z2) *
             (C(cos(mIm_z2), sin(mIm_z2))
              * w_of_z(C(y,-x))) - 1.0;
+#else
+	{
+	    cmplx zz = w_of_z(C(y,-x));
+	    double RR = exp(mRe_z2), Cv = cos(mIm_z2), Sv = sin(mIm_z2);
+	    return C(RR*(Cv*creal(zz) - Sv*cimag(zz)) - 1.0,
+		     RR*(Cv*cimag(zz) + Sv*creal(zz)));
+	}
+#endif
     }
 
     // Use Taylor series for small |z|, to avoid cancellation inaccuracy
     //   erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)
 taylor:
     {
+#ifndef __sun__
         cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
         return z * (1.1283791670955125739
                     + mz2 * (0.37612638903183752464
                              + mz2 * (0.11283791670955125739
                                       + mz2 * (0.026866170645131251760
                                                + mz2 * 0.0052239776254421878422))));
+#else
+	{
+	    double xx0, yy0, xx, yy, xxtmp;
+	    double co[5] = {
+		1.1283791670955125739,
+		0.37612638903183752464,
+		0.11283791670955125739,
+		0.026866170645131251760,
+		0.0052239776254421878422
+	    };
+	    int j;
+	    xx0 = mRe_z2; yy0 = mIm_z2; /* xx0 + I*yy0 = mz2 */
+	    xx = co[4]; yy = 0.0; /* zz = co[4] + 0i */
+	    for (j = 3; j>=0 ; j--) {
+		xxtmp = xx;
+		xx = xxtmp * xx0 - yy * yy0; 
+		yy = xxtmp * yy0 + yy * xx0;
+		xx += co[j];
+                /* zz = zz * mz2 + co[j] */
+	    }
+	    return C(xx*creal(z) - yy*cimag(z), xx*cimag(z) + yy*creal(z));
+	}
+#endif
     }
 
     /* for small |x| and small |xy|, 
@@ -290,14 +340,36 @@
     double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
     double mIm_z2 = -2*x*y; // Im(-z^2)
     if (mRe_z2 < -750) // underflow
+#ifndef __sun__
         return (x >= 0 ? 0.0 : 2.0);
+#else
+        return (x >= 0 ? C(0.0, 0.0) : C(2.0, 0.0));
+#endif
 
     if (x >= 0)
+#ifndef __sun__
         return cexp(C(mRe_z2, mIm_z2))
             * w_of_z(C(-y,x));
+#else
+    {
+	cmplx zz = w_of_z(C(-y,x));
+	double RR = exp(mRe_z2), Cv = cos(mIm_z2), Sv = sin(mIm_z2);
+	return C(RR*(Cv*creal(zz) - Sv*cimag(zz)), 
+		 RR*(Cv*cimag(zz) + Sv*creal(zz)));
+    }
+#endif
     else
+#ifndef __sun__
         return 2.0 - cexp(C(mRe_z2, mIm_z2))
             * w_of_z(C(y,-x));
+#else
+    {
+	cmplx zz = w_of_z(C(y,-x));
+	double RR = exp(mRe_z2), Cv = cos(mIm_z2), Sv = sin(mIm_z2);
+	return C(2.0 - RR*(Cv*creal(zz) - Sv*cimag(zz)), 
+		 - RR*(Cv*cimag(zz) + Sv*creal(zz)));
+    }
+#endif
 } // cerfc
 
 /******************************************************************************/
@@ -347,8 +419,18 @@
             else if (fabs(mIm_z2) < 5e-3)
                 goto taylor_realaxis;
         }
+#ifndef __sun__
         cmplx res = cexp(mz2) - w_of_z(z);
         return spi2 * C(-cimag(res), creal(res));
+#else
+	{
+	    cmplx zz = w_of_z(z);
+	    double RR = exp(creal(mz2)), Cv = cos(cimag(mz2)), 
+		Sv = sin(cimag(mz2));
+	    return C(spi2*(cimag(zz) - RR*Sv), 
+		     spi2*(RR*Cv - creal(zz)));
+	}
+#endif
     }
     else { // y < 0
         if (y > -5e-3) { // duplicate from above to avoid fabs(x) call
@@ -359,16 +441,48 @@
         }
         else if (isnan(y))
             return C(x == 0 ? 0 : NaN, NaN);
+#ifndef __sun__
         cmplx res = w_of_z(-z) - cexp(mz2);
         return spi2 * C(-cimag(res), creal(res));
+#else
+	{
+	    cmplx zz = w_of_z(C(- creal(z), - cimag(z)));
+	    double RR = exp(creal(mz2)), Cv = cos(cimag(mz2)), 
+		Sv = sin(cimag(mz2));
+	    return C(spi2*(RR*Sv - cimag(zz)), 
+		     spi2*(creal(zz) - RR*Cv));
+	}
+#endif
     }
 
     // Use Taylor series for small |z|, to avoid cancellation inaccuracy
     //     dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ...
 taylor:
+#ifndef __sun__
     return z * (1.
                 + mz2 * (0.6666666666666666666666666666666666666667
                          + mz2 * 0.2666666666666666666666666666666666666667));
+#else
+    {
+	double xx0, yy0, xx, yy, xxtmp;
+	double co[3] = {
+	    1.0,
+	    0.6666666666666666666666666666666666666667,
+	    0.2666666666666666666666666666666666666667
+	};
+	int j;
+	xx0 = creal(mz2); yy0 = cimag(mz2); /* xx0 + I*yy0 = mz2 */
+	xx = co[2]; yy = 0.0; /* zz = co[2] + 0i */
+	for (j = 1; j>=0 ; j--) {
+	    xxtmp = xx;
+	    xx = xxtmp * xx0 - yy * yy0; 
+	    yy = xxtmp * yy0 + yy * xx0;
+	    xx += co[j];
+	    /* zz = zz * mz2 + co[j] */
+	}
+	return C(xx*creal(z) - yy*cimag(z), xx*cimag(z) + yy*creal(z));
+    }
+#endif
 
     /* for small |y| and small |xy|, 
        use Taylor series to avoid cancellation inaccuracy:
@@ -417,11 +531,21 @@
                                          - 0.26666666666666666667*y2))
                          / (2*x2 - 1));
             }
+#ifndef __sun__
             return (1. / (-15 + x2*(90 + x2*(-60 + 8*x2)))) *
                 C(x * (33 + x2 * (-28 + 4*x2)
                        + y2 * (18 - 4*x2 + 4*y2)),
                   y * (-15 + x2 * (24 - 4*x2)
                        + y2 * (4*x2 - 10 - 4*y2)));
+#else
+	    {
+		double xx = 1. / (-15 + x2*(90 + x2*(-60 + 8*x2)));
+		return C(xx * (x * (33 + x2 * (-28 + 4*x2)
+				    + y2 * (18 - 4*x2 + 4*y2))),
+			 xx * (y * (-15 + x2 * (24 - 4*x2)
+				    + y2 * (4*x2 - 10 - 4*y2))));
+	    }
+#endif
         }
         else {
             double D = spi2 * im_w_of_x(x);
diff -uN libcerf-1.3/lib/w_of_z.c.ORG libcerf-1.3/lib/w_of_z.c
--- libcerf-1.3/lib/w_of_z.c.ORG	Mon Feb  4 22:04:18 2013
+++ libcerf-1.3/lib/w_of_z.c	Fri Jul 19 16:07:45 2013
@@ -79,7 +79,11 @@
 static inline cmplx cpolar(double r, double t)
 {
     if (r == 0.0 && !isnan(t))
+#ifndef __sun__
         return 0.0;
+#else
+        return C(0.0, 0.0);
+#endif
     else
         return C(r * cos(t), r * sin(t));
 }
@@ -189,7 +193,11 @@
     const double y = cimag(z);
     const double ya = fabs(y);
 
+#ifndef __sun__
     cmplx ret = 0.; // return value
+#else
+    cmplx ret = C(0.0, 0.0); // return value
+#endif
 
     double sum1 = 0, sum2 = 0, sum3 = 0, sum4 = 0, sum5 = 0;
 
@@ -264,7 +272,15 @@
             // use w(z) = 2.0*exp(-z*z) - w(-z), 
             // but be careful of overflow in exp(-z*z) 
             //                                = exp(-(xs*xs-ya*ya) -2*i*xs*ya) 
+#ifndef __sun__
             return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
+#else
+	    {
+		double RR = exp((ya-xs)*(xs+ya)), TH = 2*xs*y;
+		return C(2.0*RR*cos(TH) - creal(ret), 
+			 2.0*RR*sin(TH) - cimag(ret));
+	    }
+#endif
         }
         else
             return ret;
@@ -345,8 +361,13 @@
         if (y > 5) { // imaginary terms cancel
             faddeeva_algorithm += 10;
             const double sinxy = sin(x*y);
+#ifndef __sun__
             ret = (expx2erfcxy - c*y*sum1) * cos(2*x*y)
                 + (c*x*expx2) * sinxy * sinc(x*y, sinxy);
+#else
+            ret = C((expx2erfcxy - c*y*sum1) * cos(2*x*y)
+		    + (c*x*expx2) * sinxy * sinc(x*y, sinxy), 0.0);
+#endif
         }
         else {
             faddeeva_algorithm += 20;
@@ -368,7 +389,12 @@
         if (isnan(y))
             return C(y,y);
 
-        ret = exp(-x*x); // |y| < 1e-10, so we only need exp(-x*x) term
+#ifndef __sun__
+        ret = exp(-x*x); 
+#else
+        ret = C(exp(-x*x), 0.0); 
+#endif
+        // |y| < 1e-10, so we only need exp(-x*x) term
         // (round instead of ceil as in original paper; note that x/a > 1 here)
         double n0 = floor(x/a + 0.5); // sum in both directions, starting at n0
         double dx = a*n0 - x;
@@ -395,6 +421,11 @@
         }
     }
 finish:
+#ifndef __sun__
     return ret + C((0.5*c)*y*(sum2+sum3), 
                    (0.5*c)*copysign(sum5-sum4, creal(z)));
+#else
+    return C((0.5*c)*y*(sum2+sum3) + creal(ret), 
+	     (0.5*c)*copysign(sum5-sum4, creal(z)) + cimag(ret));
+#endif
 } // w_of_z
diff -uN libcerf-1.3/test/run_w_of_z.c.ORG libcerf-1.3/test/run_w_of_z.c
--- libcerf-1.3/test/run_w_of_z.c.ORG	Mon Feb  4 22:04:18 2013
+++ libcerf-1.3/test/run_w_of_z.c	Fri Jul 19 12:58:02 2013
@@ -38,7 +38,7 @@
     x = atof( argv[1] );
     y = atof( argv[2] );
 
-    double _Complex w = w_of_z( C(x,y) );
+    cmplx w = w_of_z( C(x,y) );
     
     double v[2][2];
     v[0][0] = creal(w);
