Merge pull request scipy#331 from pv/ticket/1741

BUG: special: replace wofz implementation with license-compatible one (scipy#1741)

scipy/special/SConscript

@@ -76,7 +76,6 @@ build_lib('cephes', '.c', 'sc_cephes')
 # F libraries
 # XXX: handle no opt flags for mach
 build_lib('mach', '.f', 'sc_mach')
-build_lib('toms', '.f', 'sc_toms')
 build_lib('amos', '.f', 'sc_amos')
 build_lib('cdflib', '.f', 'sc_cdf')
 build_lib('specfun', '.f', 'sc_specfunlib')
@@ -93,16 +92,20 @@ env.NumpyPythonExtension('orthogonal_eval', source = 'orthogonal_eval.c')
 logit_src = env.GenerateFromTemplate('_logit.c.src')
 env.NumpyPythonExtension('_logit', source = logit_src)
 
+# faddeeva extension
+env.NumpyPythonExtension('_faddeeva', source=['_faddeeva.cxx',
+                                              'faddeeva_w.cxx'])
+
 # lambertw extension
 env.NumpyPythonExtension('lambertw', source = 'lambertw.c')
 
 # Cephes extension
 src = ['_cephesmodule.c', 'amos_wrappers.c', 'specfun_wrappers.c', \
-       'toms_wrappers.c','cdf_wrappers.c','ufunc_extras.c']
+       'cdf_wrappers.c','ufunc_extras.c']
 
 env.NumpyPythonExtension('_cephes', 
                          source = src, 
-                         LIBS = ['sc_amos', 'sc_toms', 'sc_c_misc', 'sc_cephes', 'sc_mach',\
+                         LIBS = ['sc_amos', 'sc_c_misc', 'sc_cephes', 'sc_mach',\
                                  'sc_cdf', 'sc_specfunlib']) 
 
 # Specfun extension

scipy/special/__init__.py

@@ -531,6 +531,7 @@ class also has an attribute ``weights`` which returns the roots, weights,
 from spfun_stats import multigammaln
 from lambertw import lambertw
 from _logit import logit, expit
+from _faddeeva import wofz
 import add_newdocs
 
 __all__ = filter(lambda s:not s.startswith('_'),dir())

scipy/special/_cephesmodule.c

@@ -15,7 +15,6 @@
 #include "abstract.h"
 #include "cephes.h"
 #include "amos_wrappers.h"
-#include "toms_wrappers.h"
 #include "cdf_wrappers.h"
 #include "specfun_wrappers.h"
 #include "c_misc/misc.h"
@@ -206,8 +205,6 @@ static void * zetac_data[] = { (void *)zetac, (void *)zetac, };
 static void * kolmogorov_data[] = { (void *)kolmogorov, (void *)kolmogorov, };
 static void * kolmogi_data[] = { (void *)kolmogi, (void *)kolmogi, };
 
-static void * wofz_data[] = { (void *)cwofz_wrap, (void *)cwofz_wrap, };
-
 static void * besselpoly_data[] = {(void *)besselpoly, (void *)besselpoly,};
 
 static void * cdfbet3_data[] = {(void *)cdfbet3_wrap, (void *)cdfbet3_wrap};
@@ -821,10 +818,6 @@ static void Cephes_InitOperators(PyObject *dictionary) {
 	PyDict_SetItemString(dictionary, "kolmogi", f);
 	Py_DECREF(f);
 
-	f = PyUFunc_FromFuncAndData(cephes1c_functions, wofz_data, cephes_1c_types, 2, 1, 1, PyUFunc_None, "wofz", wofz_doc, 0);
-	PyDict_SetItemString(dictionary, "wofz", f);
-	Py_DECREF(f);
-
 	f = PyUFunc_FromFuncAndData(cephes3_functions, besselpoly_data, cephes_4_types, 2, 3, 1, PyUFunc_None, "besselpoly", besselpoly_doc, 0);
 	PyDict_SetItemString(dictionary, "besselpoly", f);
 	Py_DECREF(f);

scipy/special/_faddeeva.cxx

@@ -0,0 +1,108 @@
+#include <complex>
+#include <Python.h>
+#include <math.h>
+
+extern std::complex<double> Faddeeva_w(std::complex<double> z, double relerr);
+
+extern "C" {
+
+#include "numpy/npy_math.h"
+#include "numpy/ndarraytypes.h"
+#include "numpy/ufuncobject.h"
+
+#define RELERR 0    /* machine precision */
+
+static void
+wofz_loop(char **args, npy_intp *dimensions,
+          npy_intp* steps, void* data)
+{
+    npy_intp i;
+    npy_intp n = dimensions[0];
+    char *in = args[0], *out = args[1];
+    npy_intp in_step = steps[0], out_step = steps[1];
+
+    for (i = 0; i < n; i++) {
+        std::complex<double> z(((npy_cdouble*)in)->real,
+                               ((npy_cdouble*)in)->imag);
+        std::complex<double> w = Faddeeva_w(z, RELERR);
+
+        ((npy_cdouble*)out)->real = real(w);
+        ((npy_cdouble*)out)->imag = imag(w);
+
+        in += in_step;
+        out += out_step;
+    }
+}
+
+/*
+ * Definitions for the ufuncs.
+ */
+
+static PyUFuncGenericFunction wofz_funcs[1] = {&wofz_loop};
+static char types[2] = {NPY_CDOUBLE, NPY_CDOUBLE};
+static void *data[1] = {NULL};
+
+/* Module definition */
+
+static PyMethodDef module_methods[] = {
+    { NULL, NULL, 0, NULL }
+};
+
+static void _init_funcs(PyObject *m)
+{
+    PyObject *f, *d;
+
+    d = PyModule_GetDict(m);
+
+    f = PyUFunc_FromFuncAndData(wofz_funcs, data, types, 3, 1, 1,
+                                PyUFunc_None, (char*)"wofz", NULL , 0);
+    PyDict_SetItemString(d, "wofz", f);
+    Py_DECREF(f);
+}
+
+#if PY_VERSION_HEX >= 0x03000000
+
+static PyModuleDef moduledef = {
+    PyModuleDef_HEAD_INIT,
+    "_faddeeva",
+    NULL,
+    -1,
+    module_methods,
+    NULL,
+    NULL,
+    NULL,
+    NULL
+};
+
+PyMODINIT_FUNC
+PyInit__faddeeva()
+{
+    PyObject *m;
+    m = PyModule_Create(&moduledef);
+    if (!m) {
+        return NULL;
+    }
+    import_array();
+    import_umath();
+    _init_funcs(m);
+    return m;
+}
+
+#else
+
+PyMODINIT_FUNC
+init_faddeeva()
+{
+    PyObject *m;
+    m  = Py_InitModule("_faddeeva", module_methods);
+    if (m == NULL) {
+        return;
+    }
+    import_array();
+    import_umath();
+    _init_funcs(m);
+}
+
+#endif
+
+}

scipy/special/bento.info

@@ -15,9 +15,6 @@ Library:
             mach/i1mach.f,
             mach/r1mach.f,
             mach/xerror.f
-    CompiledLibrary: sc_toms
-        Sources:
-            toms/wofz.f
     CompiledLibrary: sc_amos
         Sources:
             amos/*.f
@@ -41,9 +38,11 @@ Library:
             _cephesmodule.c,
             amos_wrappers.c,
             specfun_wrappers.c,
-            toms_wrappers.c,
             cdf_wrappers.c,
             ufunc_extras.c
     Extension: specfun
          Sources:
             specfun.pyf
+    Extension: _faddeeva
+         Sources:
+            _faddeeva.cxx, faddeeva_w.cxx

scipy/special/bscript

@@ -12,3 +12,5 @@ def pre_build(context):
     context.tweak_extension("_cephes",
                           use="sc_amos sc_toms sc_c_misc sc_cephes sc_mach " \
                               "sc_cdf sc_specfunlib NPYMATH CLIB")
+
+    context.tweak_extension("_faddeeva", features="cxx cxxshlib pyext bento")