Merge pull request #8 from steppi/unholy-surgery

ENH: Add Faddeeva and SciPy special's remaining exp and log functions from

Author: steppi

include/xsf/erf.h

@@ -0,0 +1,94 @@
+#pragma once
+
+#include "faddeeva.h"
+#include "cephes/ndtr.h"
+#include "config.h"
+
+namespace xsf {
+
+inline double erf(double x) { return cephes::erf(x); }
+
+inline float erf(float x) { return erf(static_cast<double>(x)); }
+
+inline std::complex<double> erf(std::complex<double> z) { return Faddeeva::erf(z); }
+
+inline std::complex<float> erf(std::complex<float> x) {
+    return static_cast<std::complex<float>>(erf(static_cast<std::complex<double>>(x)));
+}
+
+inline double erfc(double x) { return cephes::erfc(x); }
+
+inline float erfc(float x) { return erfc(static_cast<double>(x)); }
+
+inline std::complex<double> erfc(std::complex<double> z) { return Faddeeva::erfc(z); }
+
+inline std::complex<float> erfc(std::complex<float> x) {
+    return static_cast<std::complex<float>>(erfc(static_cast<std::complex<double>>(x)));
+}
+
+inline double erfcx(double x) { return Faddeeva::erfcx(x); }
+
+inline float erfcx(float x) { return erfcx(static_cast<double>(x)); }
+
+inline std::complex<double> erfcx(std::complex<double> z) { return Faddeeva::erfcx(z); }
+
+inline std::complex<float> erfcx(std::complex<float> x) {
+    return static_cast<std::complex<float>>(erfcx(static_cast<std::complex<double>>(x)));
+}
+
+inline double erfi(double x) { return Faddeeva::erfi(x); }
+
+inline float erfi(float x) { return erfi(static_cast<double>(x)); }
+
+inline std::complex<double> erfi(std::complex<double> z) { return Faddeeva::erfi(z); }
+
+inline std::complex<float> erfi(std::complex<float> z) {
+    return static_cast<std::complex<float>>(erfi(static_cast<std::complex<double>>(z)));
+}
+
+inline double voigt_profile(double x, double sigma, double gamma) {
+    const double INV_SQRT_2 = 0.707106781186547524401;
+    const double SQRT_2PI = 2.5066282746310002416123552393401042;
+
+    if (sigma == 0) {
+        if (gamma == 0) {
+            if (std::isnan(x))
+                return x;
+            if (x == 0)
+                return INFINITY;
+            return 0;
+        }
+        return gamma / M_PI / (x * x + gamma * gamma);
+    }
+    if (gamma == 0) {
+        return 1 / SQRT_2PI / sigma * exp(-(x / sigma) * (x / sigma) / 2);
+    }
+
+    double zreal = x / sigma * INV_SQRT_2;
+    double zimag = gamma / sigma * INV_SQRT_2;
+    std::complex<double> z(zreal, zimag);
+    std::complex<double> w = Faddeeva::w(z);
+    return w.real() / sigma / SQRT_2PI;
+}
+
+inline float voigt_profile(float x, float sigma, float gamma) {
+    return voigt_profile(static_cast<double>(x), static_cast<double>(sigma), static_cast<double>(gamma));
+}
+
+inline std::complex<double> wofz(std::complex<double> z) { return Faddeeva::w(z); }
+
+inline std::complex<float> wofz(std::complex<float> x) {
+    return static_cast<std::complex<float>>(wofz(static_cast<std::complex<double>>(x)));
+}
+
+inline double dawsn(double x) { return Faddeeva::Dawson(x); }
+
+inline float dawsn(float x) { return dawsn(static_cast<double>(x)); }
+
+inline std::complex<double> dawsn(std::complex<double> z) { return Faddeeva::Dawson(z); }
+
+inline std::complex<float> dawsn(std::complex<float> x) {
+    return static_cast<std::complex<float>>(dawsn(static_cast<std::complex<double>>(x)));
+}
+
+} // namespace xsf

include/xsf/exp.h

@@ -0,0 +1,56 @@
+#pragma once
+
+#include "xsf/cephes/exp10.h"
+#include "xsf/cephes/exp2.h"
+
+namespace xsf {
+
+inline double expm1(double x) { return cephes::expm1(x); }
+
+inline float expm1(float x) { return expm1(static_cast<double>(x)); }
+
+// cexpm1(z) = cexp(z) - 1
+//
+// The imaginary part of this is easily computed via exp(z.real)*sin(z.imag)
+// The real part is difficult to compute when there is cancellation e.g. when
+// z.real = -log(cos(z.imag)).  There isn't a way around this problem  that
+// doesn't involve computing exp(z.real) and/or cos(z.imag) to higher
+// precision.
+inline std::complex<double> expm1(std::complex<double> z) {
+    if (!std::isfinite(std::real(z)) || !std::isfinite(std::imag(z))) {
+        return std::exp(z) - 1.0;
+    }
+
+    double x;
+    double ezr = 0;
+    if (std::real(z) <= -40) {
+        x = -1.0;
+    } else {
+        ezr = expm1(std::real(z));
+        x = ezr * std::cos(std::imag(z)) + cosm1(std::imag(z));
+    }
+
+    // don't compute exp(zr) too, unless necessary
+    double y;
+    if (std::real(z) > -1.0) {
+        y = (ezr + 1.0) * sin(std::imag(z));
+    } else {
+        y = exp(std::real(z)) * sin(std::imag(z));
+    }
+
+    return std::complex<double>{x, y};
+}
+
+inline std::complex<float> expm1(std::complex<float> z) {
+    return static_cast<std::complex<float>>(expm1(static_cast<std::complex<double>>(z)));
+}
+
+double exp2(double x) { return cephes::exp2(x); }
+
+float exp2(float x) { return exp2(static_cast<double>(x)); }
+
+double exp10(double x) { return cephes::exp10(x); }
+
+float exp10(float x) { return exp10(static_cast<double>(x)); }
+
+} // namespace xsf