ENH: Implement the Debye function D1(x) as a ufunc.

README.md

@@ -36,9 +36,9 @@ Most of the element-wise ufuncs are implemented by writing the core
 calculation as a templated C++ function, and using some Python code to
 automate the generation of all the necessary boilerplate and wrappers
 that implement a ufunc around the core calculation.  The exceptions
-are `logfactorial`, `log1p`, `loggamma1p`, `issnan`, and `cabssq`, which
-are implemented in C, with all boilerplate code written "by hand" in the
-C file.
+are `debye1`, `logfactorial`, `log1p`, `loggamma1p`, `issnan`, and
+`cabssq`, which are implemented in C, with all boilerplate code written
+"by hand" in the C file.
 
 | Function                                      | Description                                                   |
 | --------                                      | -----------                                                   |
@@ -56,6 +56,7 @@ C file.
 | [`smoothstep5`](#smoothstep5)                 | Smooth step using a degree 5 polynomial                       |
 | [`trapezoid_pulse`](#trapezoid_pulse)         | Trapezoid pulse function                                      |
 | [`pow1pm1`](#pow1pm1)                         | Compute `(1 + x)**y - 1`                                      |
+| [`debye1`](#debye1)                           | Compute the Debye function D1(x)                              |
 | [`expint1`](#expint1)                         | Exponential integral E₁ for real inputs                       |
 | [`log1p_doubledouble`](#log1p_doubledouble)   | `log(1 + z)` for complex z.                                   |
 | [`log1p_theorem4`](#log1p_theorem4)           | `log(1 + z)` for complex z.                                   |
@@ -328,7 +329,7 @@ the script `trapezoid_pulse_demo.py` in the `examples` directory):
 
 ![trapezoid_pulse plot1](https://github.com/WarrenWeckesser/ufunclab/blob/main/examples/trapezoid_pulse_demo.png)
 
-### `pow1pm1`
+#### `pow1pm1`
 
 `pow1pm1(x, y)` computes `(1 + x)**y - 1` for `x >= -1`.
 
@@ -352,6 +353,25 @@ The naive calculation provides less than six digits of precision:
 
 ```
 
+#### `debye1`
+
+`debye1(x)` computes the Debye function D1(x).
+
+See the wikipedia article
+
+    https://en.wikipedia.org/wiki/Debye_function
+
+for more details.
+
+```
+>>> from ufunclab import debye1
+
+>>> debye1([-2, -1.5, 0, 0.25, 0.5, 2.5, 50, 100])
+array([1.60694728, 1.43614531, 1.        , 0.93923503, 0.88192716,
+       0.53878957, 0.03289868, 0.01644934])
+
+```
+
 #### `expint1`
 
 `expint1(x)` computes the exponential integral E₁ for the real input x.

src/debye1/check_debye1_generated.c

@@ -0,0 +1,27 @@
+//
+// A simple C main program to exercise the code in debye1_generated.c.
+// Compile with (for example):
+//     $ gcc check_debye1_generated.c debye1_generated.c -o check_debye1_generated
+// Pass x values on the command line, e.g. on Linux:
+//     $ ./check_debye1_generated 0 0.25 5.5 -12.0
+// That will print:
+//     0.00000000000000000e+00  1.00000000000000000e+00
+//     2.50000000000000000e-01  9.39235027193614513e-01
+//     5.50000000000000000e+00  2.94239966231542471e-01
+//    -1.20000000000000000e+01  6.13707118265430740e+00
+//
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "debye1_generated.h"
+
+
+int main(int argc, char *argv[])
+{
+    for (int i = 1; i < argc; ++i) {
+        double x = strtod(argv[i], NULL);
+        double y = debye1(x);
+        printf("%25.17e  %.17e\n", x, y);
+    }
+}

src/debye1/debye1_generated.h

@@ -0,0 +1,10 @@
+
+// Do not edit this file!
+// This file was generated by a script.
+
+#ifndef DEBYE1_GENERATED_H
+#define DEBYE1_GENERATED_H
+
+double debye1(const double);
+
+#endif

src/debye1/debye1_ufunc.c

@@ -0,0 +1,139 @@
+//
+// debye1_ufunc.c
+//
+
+#define PY_SSIZE_T_CLEAN
+#include "Python.h"
+
+#define NPY_NO_DEPRECATED_API NPY_API_VERSION
+#include "numpy/ndarraytypes.h"
+#include "numpy/ufuncobject.h"
+
+#include "debye1_generated.h"
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// The ufunc "inner loops".
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+static void debye1_double_loop(char **args, const npy_intp *dimensions,
+                               const npy_intp* steps, void* data)
+{
+    char *in = args[0];
+    char *out = args[1];
+    npy_intp in_step = steps[0];
+    npy_intp out_step = steps[1];
+
+    for (npy_intp i = 0; i < dimensions[0]; ++i, in += in_step, out += out_step) {
+        double x = *(double *)in;
+        *((double *)out) = debye1(x);
+    }
+}
+
+static void debye1_float_loop(char **args, const npy_intp *dimensions,
+                                     const npy_intp* steps, void* data)
+{
+    char *in = args[0];
+    char *out = args[1];
+    npy_intp in_step = steps[0];
+    npy_intp out_step = steps[1];
+
+    for (npy_intp i = 0; i < dimensions[0]; ++i, in += in_step, out += out_step) {
+        double x = (double) *(float *)in;
+        *((float *)out) = debye1(x);
+    }
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// ufunc configuration data.
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+// This array of loop function pointers will be passed to PyUFunc_FromFuncAndData,
+// along with the arrays debye1_typecodes and debye1_data.
+PyUFuncGenericFunction debye1_funcs[] = {
+    (PyUFuncGenericFunction) &debye1_float_loop,
+    (PyUFuncGenericFunction) &debye1_double_loop
+};
+
+#define DEBYE1_NLOOPS \
+    (sizeof(debye1_funcs) / sizeof(debye1_funcs[0]))
+
+// These are the input and return type codes for the two loops.  It is
+// created as a 1-d C array, but it can be interpreted as a 2-d array with
+// shape (num_loops, num_args).  (num_args is the sum of the number of
+// input args and output args.  In this case num_args is 2.)
+static char debye1_typecodes[] = {
+    NPY_FLOAT, NPY_FLOAT,
+    NPY_DOUBLE, NPY_DOUBLE
+};
+
+// The ufunc will not use the 'data' array.
+static void *debye1_data[] = {NULL, NULL};
+
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Python extension module definitions.
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+
+static PyMethodDef Debye1Methods[] = {
+        {NULL, NULL, 0, NULL}
+};
+
+static struct PyModuleDef moduledef = {
+    PyModuleDef_HEAD_INIT,
+    .m_name = "_debye1",
+    .m_doc = "Module that defines the debye1 function.",
+    .m_size = -1,
+    .m_methods = Debye1Methods
+};
+
+
+#define DEBYE1_DOCSTRING \
+"debye1(x, /, ...)\n"                                 \
+"\n"                                                  \
+"The Debye function D1(x).\n"                         \
+"\n"                                                  \
+"See https://en.wikipedia.org/wiki/Debye_function\n"  \
+"\n"
+
+
+PyMODINIT_FUNC PyInit__debye1(void)
+{
+    PyObject *module;
+    PyObject *debye1_ufunc;
+    int nin, nout;
+    int status;
+
+    module = PyModule_Create(&moduledef);
+    if (!module) {
+        return NULL;
+    }
+
+    import_array();
+    import_umath();
+
+    // Create the debye1 ufunc object.
+    nin = 1;
+    nout = 1;
+    debye1_ufunc = PyUFunc_FromFuncAndData(debye1_funcs,
+                                           debye1_data,
+                                           debye1_typecodes,
+                                           DEBYE1_NLOOPS, nin, nout,
+                                           PyUFunc_None,
+                                           "debye1",
+                                           DEBYE1_DOCSTRING, 0);
+    if (debye1_ufunc == NULL) {
+        Py_DECREF(module);
+        return NULL;
+    }
+
+    // Add the ufunc to the module.
+    status = PyModule_AddObject(module, "debye1",
+                                (PyObject *) debye1_ufunc);
+    if (status == -1) {
+        Py_DECREF(debye1_ufunc);
+        Py_DECREF(module);
+        return NULL;
+    }
+
+    return module;
+}

src/debye1/generate_debye1_c.py

@@ -0,0 +1,134 @@
+"""
+This script generates the files debye1_generated.c and debye1_generated.h.
+"""
+
+import mpmath
+from mpmath import mp
+
+
+def mp_debye1(x):
+    """
+    Compute the Debye function D1(x).
+
+    mp.quad is used, so it might be necessary to set mp.dps to a value
+    much larger than the desired output precision.  For example, mp.dps=150
+    is required to get an accurate double precision result from
+    mp_debye1(1e100).
+    """
+    if x == 0:
+        return mp.one
+
+    def integrand(t):
+        if t == 0:
+            return mp.one
+        return t / mp.expm1(t)
+
+    return mp.quad(integrand, [0, x])/x
+
+
+_preamble = """
+// Do not edit this file!
+// This file was generated by a script.
+"""
+
+
+def generate_polynomial_function(f, name, coeffs, scale_by_x=False,
+                                 comment=None):
+    print(f"static double {name}(const double x)", file=f)
+    print("{", file=f)
+    if comment is not None:
+        lines = comment.splitlines()
+        for line in lines:
+            print(f'    // {line}', file=f)
+    print("    const double y = ", end='', file=f)
+    nterms = len(coeffs)
+    if scale_by_x:
+        print("(", end='', file=f)
+    for i, mp_coeff in enumerate(coeffs[::-1]):
+        coeff = float(mp_coeff)
+        print(coeff, end='', file=f)
+        if i == nterms - 1:
+            print(")"*(nterms - 1), end='', file=f)
+            if scale_by_x:
+                print(")/x", end='', file=f)
+            print(";", file=f)
+        else:
+            print(" + x*(", end='', file=f)
+    print("    return y;", file=f)
+    print("}", file=f)
+
+
+def generate_code(f, intervals):
+    mp.dps = 200
+
+    print(_preamble, file=f)
+    print(file=f)
+    print('#define DEBYE1_ASYMP_CONST 1.6449340668482264  // pi**2/6', file=f)
+    print(file=f)
+    for i, (a, b, n, order, scale_by_x) in enumerate(intervals):
+        edges = mpmath.linspace(a, b, n + 1)
+        for k in range(n):
+            x0 = edges[k]
+            x1 = edges[k+1]
+            if scale_by_x:
+                poly, err = mp.chebyfit(lambda t: t*mp_debye1(t), [x0, x1],
+                                        order, error=True)
+                print(f"x0={float(x0):8.4f} x1={float(x1):8.4f} "
+                      f"err={float(err/x0):8.3e}")
+            else:
+                poly, err = mp.chebyfit(lambda t: mp_debye1(t), [x0, x1],
+                                        order, error=True)
+                print(f"x0={float(x0):8.4f} x1={float(x1):8.4f} "
+                      f"err={float(err):8.3e}")
+            funcname = f'debye1_{i}_{k:02}'
+            comment = f'Approximate debye1 on the interval [{x0}, {x1}].'
+            generate_polynomial_function(f, funcname, poly, scale_by_x,
+                                         comment=comment)
+            print(file=f)
+
+    print('double debye1(const double x)', file=f)
+    print('{', file=f)
+    print('    if (x < 0) {', file=f)
+    print('        return debye1(-x) - x/2;', file=f)
+    print('    }', file=f)
+    for i, (a, b, n, order, scale_by_x) in enumerate(intervals):
+        print(f'    if (x < {b}) {{', file=f)
+        print(f'        int k = (int) ({n}*((x - {a})/{b - a}));', file=f)
+        print('        switch (k) {', file=f)
+        for k in range(n):
+            print(f'            case {k:2}: return debye1_{i}_{k:02}(x);',
+                  file=f)
+        print('        }', file=f)
+        print('    }', file=f)
+    print(f'    // x >= {b}', file=f)
+    print('    return DEBYE1_ASYMP_CONST/x;', file=f)
+    print('}', file=f)
+
+
+if __name__ == "__main__":
+
+    # Each interval in this list will be subdivided into n subintervals,
+    # and in each subinterval, an approximation will be generated by
+    # mp.chebyfit with the given order.  If scale_by_x is True, chebyfit
+    # is applied to mp_debye1(x)*x.
+    intervals = [
+        #  x0, x1, n, order, scale_by_x
+        ( 0.0, 10.0, 50, 8, False),
+        (10.0, 20.0, 50, 8, True),
+        (20.0, 30.0, 10, 7, True),
+        (30.0, 40.0, 10, 4, True),
+    ]
+
+    fnamebase = 'debye1_generated'
+    print(f'Generating {fnamebase}.c')
+    with open(fnamebase + '.c', 'w') as f:
+        generate_code(f, intervals)
+    print(f'Generating {fnamebase}.h')
+    with open(fnamebase + '.h', 'w') as f:
+        print(_preamble, file=f)
+        print('#ifndef DEBYE1_GENERATED_H', file=f)
+        print('#define DEBYE1_GENERATED_H', file=f)
+        print(file=f)
+        print('double debye1(const double);', file=f)
+        print(file=f)
+        print('#endif', file=f)

ufunclab/__init__.py

@@ -8,6 +8,7 @@
 from ._abs_squared import abs_squared
 from ._cabssq import cabssq
 from ._log1p import log1p_theorem4, log1p_doubledouble
+from ._debye1 import debye1
 from ._expint1 import expint1, logexpint1
 from ._pow1pm1 import pow1pm1
 from ._logistic import logistic, logistic_deriv, log_logistic, swish

ufunclab/meson.build

@@ -269,6 +269,10 @@ endforeach
 #----------------------------------------------------------------------
 
 modules = [
+    ['_debye1',
+        ['../src/debye1/debye1_ufunc.c',
+         '../src/debye1/debye1_generated.h',
+         '../src/debye1/debye1_generated.c']],
     ['_cabssq',
         ['../src/cabssq/cabssq_ufunc.c']],
     ['_log1p',