add complex error functions (via http://ab-initio.mit.edu/Faddeeva)

Makefile

@@ -42,8 +42,8 @@ JL_LIBS = julia-release julia-debug
 # private libraries, that are installed in $(PREFIX)/lib/julia
 JL_PRIVATE_LIBS = amd arpack cholmod colamd fftw3 fftw3f fftw3_threads \
                   fftw3f_threads glpk glpk_wrapper gmp gmp_wrapper grisu \
-                  history openlibm pcre random readline Rmath spqr \
-                  suitesparse_wrapper tk_wrapper umfpack z openblas
+                  history Faddeeva_wrapper openlibm pcre random readline \
+	          Rmath spqr suitesparse_wrapper tk_wrapper umfpack z openblas
 
 PREFIX ?= julia-$(JULIA_COMMIT)
 install: release

base/export.jl

@@ -393,13 +393,16 @@ export
     csc,
     cscd,
     csch,
+    dawson,
     degrees2radians,
     den,
     digamma,
     div,
     eps,
     erf,
     erfc,
+    erfcx,
+    erfi,
     exp,
     exp2,
     expm1,

base/math.jl

@@ -7,7 +7,8 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,
        acosd, acotd, acscd, asecd, asind, atand, atan2,
        radians2degrees, degrees2radians,
        log, log2, log10, log1p, logb, exp, exp2, expm1, 
-       cbrt, sqrt, square, erf, erfc, ceil, floor, trunc, round, 
+       cbrt, sqrt, square, erf, erfc, erfcx, erfi, dawson,
+       ceil, floor, trunc, round, 
        lgamma, hypot, gamma, lfact, max, min, ilogb, ldexp, frexp,
        clamp, modf, ^, 
        airy, airyai, airyprime, airyaiprime, airybi, airybiprime,
@@ -623,4 +624,34 @@ function zeta(z::Number)
     eta(z) * zz/(zz-2)
 end
 
+# wrappers for complex Faddeeva functions; these will get a lot simpler,
+# and can call openlibm directly, once ccall supports C99 complex types.
+let
+    for f in (:erf, :erfc, :erfcx, :erfi, :Dawson)
+        @eval begin
+            global $f
+            function ($f)(z::Complex128)
+                local w::Array{Float64,1} = Array(Float64,2)
+                ccall(($(string("wrapFaddeeva_",f)),:libFaddeeva_wrapper), Void, (Ptr{Complex128},Ptr{Complex128},Float64,), pointer(w), &z, zero(Float64))
+                return complex128(w[1],w[2])
+            end
+            function ($f)(z::Complex64)
+                local w::Array{Float64,1} = Array(Float64,2)
+                ccall(($(string("wrapFaddeeva_",f)),:libFaddeeva_wrapper), Void, (Ptr{Complex128},Ptr{Complex128},Float64,), pointer(w), &complex128(z), float64(eps(Float32)))
+                return complex64(w[1],w[2])
+            end
+            ($f)(z::Complex) = ($f)(complex128(z))
+        end
+    end
+end
+for f in (:erfcx, :erfi, :Dawson)
+    @eval begin
+        ($f)(x::Float64) = ccall(($(string("Faddeeva_",f,"_re")),:libopenlibm), Float64, (Float64,), x)
+        ($f)(x::Float32) = float32(ccall(($(string("Faddeeva_",f,"_re")),:libopenlibm), Float64, (Float64,), float64(x)))
+        ($f)(x::Real) = ($f)(float(x))
+        @vectorize_1arg Number $f
+    end
+end
+dawson(z) = Dawson(z) # convert to lower-case naming convention
+
 end # module

deps/Faddeeva_wrapper.c

@@ -0,0 +1,45 @@
+/* Copyright (c) 2012 Massachusetts Institute of Technology
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ * 
+ * The above copyright notice and this permission notice shall be
+ * included in all copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 
+ */
+
+/* C89-compatible wrappers for complex Faddeeva functions, that 
+   pass and return complex values via a double* which points to
+   the real and imaginary parts consecutively.  (Note that this
+   is binary-compatible with both C99 complex numbers and with
+   C++ std::complex<double>.) */
+
+#include "Faddeeva.h"
+
+#define WRAP(func)							\
+void wrap ## Faddeeva_ ## func(double *func, const double *z, double relerr) \
+{									\
+     const double complex *zc = (const double complex *) z;		\
+     double complex c = Faddeeva_ ## func(*zc, relerr);			\
+     func[0] = creal(c);						\
+     func[1] = cimag(c);						\
+}
+
+WRAP(w)
+WRAP(erfcx)
+WRAP(erf)
+WRAP(erfi)
+WRAP(erfc)
+WRAP(Dawson)

deps/Makefile

@@ -13,7 +13,7 @@ include $(JULIAHOME)/Make.inc
 unexport CONFIG_SITE
 
 STAGE1_DEPS = uv openlibm random rmath double-conversion glpk
-STAGE2_DEPS = gmp-wrapper
+STAGE2_DEPS = gmp-wrapper Faddeeva-wrapper
 STAGE3_DEPS = suitesparse-wrapper
 
 ifeq ($(OS), Linux)
@@ -384,13 +384,24 @@ $(OPENLIBM_OBJ_SOURCE): openlibm/Makefile
 	$(MAKE) -C openlibm $(OPENLIBM_FLAGS) CC="$(CC)" FC="$(FC)" USECLANG=$(USECLANG) USEGCC=$(USEGCC)
 $(OPENLIBM_OBJ_TARGET): $(OPENLIBM_OBJ_SOURCE) | $(BUILD)/lib
 	cp $< $@
-install-openlibm: $(OPENLIBM_OBJ_TARGET)
+install-openlibm: $(OPENLIBM_OBJ_TARGET) $(BUILD)/lib/libFaddeeva_wrapper.$(SHLIB_EXT)
 
 clean-openlibm:
 	-$(MAKE) -C openlibm distclean
 	-rm $(OPENLIBM_OBJ_TARGET)
 distclean-openlibm: clean-openlibm      
 
+# Wrapper for openlibm/Faddeeva since ccall doesn't support C99 complex:
+$(BUILD)/lib/libFaddeeva_wrapper.$(SHLIB_EXT): Faddeeva_wrapper.c $(OPENLIBM_OBJ_TARGET)
+	$(CC) $(CPPFLAGS) $(CFLAGS) $(LDFLAGS) -O2 -shared $(fPIC) -I openlibm/Faddeeva Faddeeva_wrapper.c -o $@ -L $(BUILD)/lib -lopenlibm
+	$(INSTALL_NAME_CMD)libFaddeeva_wrapper.$(SHLIB_EXT) $@
+	touch $@
+install-Faddeeva-wrapper: $(BUILD)/lib/libFaddeeva_wrapper.$(SHLIB_EXT)
+
+clean-Faddeeva-wrapper:
+	-rm -f $(OPENLIBM_OBJ_TARGET) $(BUILD)/lib/libFaddeeva_wrapper.$(SHLIB_EXT)
+distclean-Faddeeva-wrapper: clean-Faddeeva-wrapper
+
 ## Rmath ##
 
 RMATH_OBJ_TARGET = $(BUILD)/lib/libRmath.$(SHLIB_EXT)

doc/stdlib/base.rst

@@ -1030,6 +1030,21 @@ Mathematical Functions
 
    Compute the complementary error function of ``x``
 
+.. function:: erfcx(x)
+
+   Compute the scaled complementary error function of ``x``,
+   defined by :math:`e^{x^2} \operatorname{erfc}(x)`
+
+.. function:: erfi(x)
+
+   Compute the imaginary error function of ``x``,
+   defined by :math:`-i \operatorname{erf}(ix)`
+
+.. function:: dawson(x)
+
+   Compute the Dawson function (scaled imaginary error function) of ``x``,
+   defined by :math:`\frac{\pi}{2} e^{-x^2} \operatorname{erfi}(x)`
+
 .. function:: real(z)
 
    Return the real part of the complex number ``z``

test/math.jl

@@ -1,3 +1,15 @@
+# error functions
+@assert_approx_eq erf(1) 0.84270079294971486934
+@assert_approx_eq erf(1+2im) -0.53664356577856503399-5.0491437034470346695im
+@assert_approx_eq erfc(1) 0.15729920705028513066
+@assert_approx_eq erfc(1+2im) 1.5366435657785650340+5.0491437034470346695im
+@assert_approx_eq erfcx(1) 0.42758357615580700442
+@assert_approx_eq erfcx(1+2im) 0.14023958136627794370-0.22221344017989910261im
+@assert_approx_eq erfi(1) 1.6504257587975428760
+@assert_approx_eq erfi(1+2im) -0.011259006028815025076+1.0036063427256517509im
+@assert_approx_eq dawson(1) 0.53807950691276841914
+@assert_approx_eq dawson(1+2im) -13.388927316482919244-11.828715103889593303im
+
 # airy
 @assert_approx_eq airy(1.8) 0.0470362168668458052247
 @assert_approx_eq airyprime(1.8) -0.0685247801186109345638