summaryrefslogtreecommitdiff
path: root/source/luametatex/source/libraries/libcerf/err_fcts.c
diff options
context:
space:
mode:
Diffstat (limited to 'source/luametatex/source/libraries/libcerf/err_fcts.c')
-rw-r--r--source/luametatex/source/libraries/libcerf/err_fcts.c438
1 files changed, 438 insertions, 0 deletions
diff --git a/source/luametatex/source/libraries/libcerf/err_fcts.c b/source/luametatex/source/libraries/libcerf/err_fcts.c
new file mode 100644
index 000000000..9c0c7aed9
--- /dev/null
+++ b/source/luametatex/source/libraries/libcerf/err_fcts.c
@@ -0,0 +1,438 @@
+/* Library libcerf:
+ * Compute complex error functions, based on a new implementation of
+ * Faddeeva's w_of_z. Also provide Dawson and Voigt functions.
+ *
+ * File err_fcts.c:
+ * Computate Dawson, Voigt, and several error functions,
+ * based on erfcx, im_w_of_x, w_of_z as implemented in separate files.
+ *
+ * Given w(z), the error functions are mostly straightforward
+ * to compute, except for certain regions where we have to
+ * switch to Taylor expansions to avoid cancellation errors
+ * [e.g. near the origin for erf(z)].
+ *
+ * Copyright:
+ * (C) 2012 Massachusetts Institute of Technology
+ * (C) 2013 Forschungszentrum Jülich GmbH
+ *
+ * Licence:
+ * 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.
+ *
+ * Authors:
+ * Steven G. Johnson, Massachusetts Institute of Technology, 2012, core author
+ * Joachim Wuttke, Forschungszentrum Jülich, 2013, package maintainer
+ *
+ * Website:
+ * http://apps.jcns.fz-juelich.de/libcerf
+ *
+ * Revision history:
+ * ../CHANGELOG
+ *
+ * Man pages:
+ * cerf(3), dawson(3), voigt(3)
+ */
+
+#include "cerf.h"
+#include <math.h>
+#include "defs.h" // defines _cerf_cmplx, NaN, C, cexp, ...
+
+const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2
+const double s2pi = 2.5066282746310005024157652848110; // sqrt(2*pi)
+const double pi = 3.141592653589793238462643383279503;
+
+/******************************************************************************/
+/* Simple wrappers: cerfcx, cerfi, erfi, dawson */
+/******************************************************************************/
+
+_cerf_cmplx cerfcx(_cerf_cmplx z)
+{
+ // Compute erfcx(z) = exp(z^2) erfc(z),
+ // the complex underflow-compensated complementary error function,
+ // trivially related to Faddeeva's w_of_z.
+
+ return w_of_z(C(-cimag(z), creal(z)));
+}
+
+_cerf_cmplx cerfi(_cerf_cmplx z)
+{
+ // Compute erfi(z) = -i erf(iz),
+ // the rotated complex error function.
+
+ _cerf_cmplx e = cerf(C(-cimag(z),creal(z)));
+ return C(cimag(e), -creal(e));
+}
+
+double erfi(double x)
+{
+ // Compute erfi(x) = -i erf(ix),
+ // the imaginary error function.
+
+ return x*x > 720 ? (x > 0 ? Inf : -Inf) : exp(x*x) * im_w_of_x(x);
+}
+
+double dawson(double x)
+{
+
+ // Compute dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x),
+ // Dawson's integral for a real argument.
+
+ return spi2 * im_w_of_x(x);
+}
+
+double re_w_of_z( double x, double y )
+{
+ return creal( w_of_z( C(x,y) ) );
+}
+
+double im_w_of_z( double x, double y )
+{
+ return cimag( w_of_z( C(x,y) ) );
+}
+
+/******************************************************************************/
+/* voigt */
+/******************************************************************************/
+
+double voigt( double x, double sigma, double gamma )
+{
+ // Joachim Wuttke, January 2013.
+
+ // Compute Voigt's convolution of a Gaussian
+ // G(x,sigma) = 1/sqrt(2*pi)/|sigma| * exp(-x^2/2/sigma^2)
+ // and a Lorentzian
+ // L(x,gamma) = |gamma| / pi / ( x^2 + gamma^2 ),
+ // namely
+ // voigt(x,sigma,gamma) =
+ // \int_{-infty}^{infty} dx' G(x',sigma) L(x-x',gamma)
+ // using the relation
+ // voigt(x,sigma,gamma) = Re{ w(z) } / sqrt(2*pi) / |sigma|
+ // with
+ // z = (x+i*|gamma|) / sqrt(2) / |sigma|.
+
+ // Reference: Abramowitz&Stegun (1964), formula (7.4.13).
+
+ double gam = gamma < 0 ? -gamma : gamma;
+ double sig = sigma < 0 ? -sigma : sigma;
+
+ if ( gam==0 ) {
+ if ( sig==0 ) {
+ // It's kind of a delta function
+ return x ? 0 : Inf;
+ } else {
+ // It's a pure Gaussian
+ return exp( -x*x/2/(sig*sig) ) / s2pi / sig;
+ }
+ } else {
+ if ( sig==0 ) {
+ // It's a pure Lorentzian
+ return gam / pi / (x*x + gam*gam);
+ } else {
+ // Regular case, both parameters are nonzero
+ _cerf_cmplx z = complex_mul_cr(C(x, gam), 1. / sqrt(2) / sig);
+ return creal( w_of_z(z) ) / s2pi / sig;
+ // TODO: correct and activate the following:
+// double w = sqrt(gam*gam+sig*sig); // to work in reduced units
+// _cerf_cmplx z = C(x/w,gam/w) / sqrt(2) / (sig/w);
+// return creal( w_of_z(z) ) / s2pi / (sig/w);
+ }
+ }
+}
+
+/******************************************************************************/
+/* cerf */
+/******************************************************************************/
+
+_cerf_cmplx cerf(_cerf_cmplx z)
+{
+
+ // Steven G. Johnson, October 2012.
+
+ // Compute erf(z), the complex error function,
+ // using w_of_z except for certain regions.
+
+ double x = creal(z), y = cimag(z);
+
+ if (y == 0)
+ return C(erf(x), y); // preserve sign of 0
+ if (x == 0) // handle separately for speed & handling of y = Inf or NaN
+ return C(x, // preserve sign of 0
+ /* handle y -> Inf limit manually, since
+ exp(y^2) -> Inf but Im[w(y)] -> 0, so
+ IEEE will give us a NaN when it should be Inf */
+ y*y > 720 ? (y > 0 ? Inf : -Inf)
+ : exp(y*y) * im_w_of_x(y));
+
+ 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
+ return (x >= 0 ? C(1.0, 0.0) : C(-1.0, 0.0));;
+
+ /* Handle positive and negative x via different formulas,
+ using the mirror symmetries of w, to avoid overflow/underflow
+ problems from multiplying exponentially large and small quantities. */
+ if (x >= 0) {
+ if (x < 8e-2) {
+ if (fabs(y) < 1e-2)
+ goto taylor;
+ else if (fabs(mIm_z2) < 5e-3 && x < 5e-3)
+ goto taylor_erfi;
+ }
+ /* don't use complex exp function, since that will produce spurious NaN
+ values when multiplying w in an overflow situation. */
+ return complex_sub_rc(1.0, complex_mul_rc(exp(mRe_z2), complex_mul_cc(C(cos(mIm_z2), sin(mIm_z2)), w_of_z(C(-y, x)))));
+ }
+ else { // x < 0
+ if (x > -8e-2) { // duplicate from above to avoid fabs(x) call
+ if (fabs(y) < 1e-2)
+ goto taylor;
+ else if (fabs(mIm_z2) < 5e-3 && x > -5e-3)
+ goto taylor_erfi;
+ }
+ else if (isnan(x))
+ 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. */
+ return complex_add_rc(-1.0, complex_mul_rc(exp(mRe_z2), complex_mul_cc(C(cos(mIm_z2), sin(mIm_z2)), w_of_z(C(y, -x)))));
+
+ }
+
+ // 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:
+ {
+ _cerf_cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
+ return
+ complex_mul_cc(z, complex_add_rc(1.1283791670955125739,
+ complex_mul_cc(mz2, complex_add_rc(0.37612638903183752464,
+ complex_mul_cc(mz2, complex_add_rc(0.11283791670955125739,
+ complex_mul_cc(mz2, complex_add_rc(0.026866170645131251760,
+ complex_mul_cr(mz2, 0.0052239776254421878422)))))))));
+
+
+ }
+
+ /* for small |x| and small |xy|,
+ use Taylor series to avoid cancellation inaccuracy:
+ erf(x+iy) = erf(iy)
+ + 2*exp(y^2)/sqrt(pi) *
+ [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ...
+ - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
+ where:
+ erf(iy) = exp(y^2) * Im[w(y)]
+ */
+taylor_erfi:
+ {
+ double x2 = x*x, y2 = y*y;
+ double expy2 = exp(y2);
+ return C
+ (expy2 * x * (1.1283791670955125739
+ - x2 * (0.37612638903183752464
+ + 0.75225277806367504925*y2)
+ + x2*x2 * (0.11283791670955125739
+ + y2 * (0.45135166683820502956
+ + 0.15045055561273500986*y2))),
+ expy2 * (im_w_of_x(y)
+ - x2*y * (1.1283791670955125739
+ - x2 * (0.56418958354775628695
+ + 0.37612638903183752464*y2))));
+ }
+} // cerf
+
+/******************************************************************************/
+/* cerfc */
+/******************************************************************************/
+
+_cerf_cmplx cerfc(_cerf_cmplx z)
+{
+ // Steven G. Johnson, October 2012.
+
+ // Compute erfc(z) = 1 - erf(z), the complex complementary error function,
+ // using w_of_z except for certain regions.
+
+ double x = creal(z), y = cimag(z);
+
+ if (x == 0.)
+ return C(1,
+ /* handle y -> Inf limit manually, since
+ exp(y^2) -> Inf but Im[w(y)] -> 0, so
+ IEEE will give us a NaN when it should be Inf */
+ y*y > 720 ? (y > 0 ? -Inf : Inf)
+ : -exp(y*y) * im_w_of_x(y));
+ if (y == 0.) {
+ if (x*x > 750) // underflow
+ return C(x >= 0 ? 0.0 : 2.0,
+ -y); // preserve sign of 0
+ return C(x >= 0 ? exp(-x*x) * erfcx(x)
+ : 2. - exp(-x*x) * erfcx(-x),
+ -y); // preserve sign of zero
+ }
+
+ 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
+ return C((x >= 0 ? 0.0 : 2.0), 0.0);
+
+ if (x >= 0)
+ return cexp(complex_mul_cc(C(mRe_z2, mIm_z2), w_of_z(C(-y,x))));
+ else
+ return complex_sub_rc(2.0, complex_mul_cc(cexp(C(mRe_z2, mIm_z2)), w_of_z(C(y, -x))));
+} // cerfc
+
+/******************************************************************************/
+/* cdawson */
+/******************************************************************************/
+
+_cerf_cmplx cdawson(_cerf_cmplx z)
+{
+
+ // Steven G. Johnson, October 2012.
+
+ // Compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z),
+ // Dawson's integral for a complex argument,
+ // using w_of_z except for certain regions.
+
+ double x = creal(z), y = cimag(z);
+
+ // handle axes separately for speed & proper handling of x or y = Inf or NaN
+ if (y == 0)
+ return C(spi2 * im_w_of_x(x),
+ -y); // preserve sign of 0
+ if (x == 0) {
+ double y2 = y*y;
+ if (y2 < 2.5e-5) { // Taylor expansion
+ return C(x, // preserve sign of 0
+ y * (1.
+ + y2 * (0.6666666666666666666666666666666666666667
+ + y2 * 0.26666666666666666666666666666666666667)));
+ }
+ return C(x, // preserve sign of 0
+ spi2 * (y >= 0
+ ? exp(y2) - erfcx(y)
+ : erfcx(-y) - exp(y2)));
+ }
+
+ double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
+ double mIm_z2 = -2*x*y; // Im(-z^2)
+ _cerf_cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
+
+ /* Handle positive and negative x via different formulas,
+ using the mirror symmetries of w, to avoid overflow/underflow
+ problems from multiplying exponentially large and small quantities. */
+ if (y >= 0) {
+ if (y < 5e-3) {
+ if (fabs(x) < 5e-3)
+ goto taylor;
+ else if (fabs(mIm_z2) < 5e-3)
+ goto taylor_realaxis;
+ }
+ _cerf_cmplx res = complex_sub_cc(cexp(mz2), w_of_z(z));
+ return complex_mul_rc(spi2, C(-cimag(res), creal(res)));
+ }
+ else { // y < 0
+ if (y > -5e-3) { // duplicate from above to avoid fabs(x) call
+ if (fabs(x) < 5e-3)
+ goto taylor;
+ else if (fabs(mIm_z2) < 5e-3)
+ goto taylor_realaxis;
+ }
+ else if (isnan(y))
+ return C(x == 0 ? 0 : NaN, NaN);
+ {
+ _cerf_cmplx res = complex_sub_cc(w_of_z(complex_neg(z)), cexp(mz2));
+ return complex_mul_rc(spi2, C(-cimag(res), creal(res)));
+ }
+ }
+
+ // Use Taylor series for small |z|, to avoid cancellation inaccuracy
+ // dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ...
+taylor:
+ return complex_mul_cc(z, complex_add_rc(1.,
+ complex_mul_cc(mz2, complex_add_rc(0.6666666666666666666666666666666666666667,
+ complex_mul_cr(mz2, 0.2666666666666666666666666666666666666667)))));
+ /* for small |y| and small |xy|,
+ use Taylor series to avoid cancellation inaccuracy:
+ dawson(x + iy)
+ = D + y^2 (D + x - 2Dx^2)
+ + y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3)
+ + iy [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3)
+ + y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ...
+ where D = dawson(x)
+
+ However, for large |x|, 2Dx -> 1 which gives cancellation problems in
+ this series (many of the leading terms cancel). So, for large |x|,
+ we need to substitute a continued-fraction expansion for D.
+
+ dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...))))))
+
+ The 6 terms shown here seems to be the minimum needed to be
+ accurate as soon as the simpler Taylor expansion above starts
+ breaking down. Using this 6-term expansion, factoring out the
+ denominator, and simplifying with Maple, we obtain:
+
+ Re dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / x
+ = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4
+ Im dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / y
+ = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4
+
+ Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction
+ expansion for the real part, and a 2-term expansion for the imaginary
+ part. (This avoids overflow problems for huge |x|.) This yields:
+
+ Re dawson(x + iy) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x)
+ Im dawson(x + iy) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1)
+
+ */
+taylor_realaxis:
+ {
+ double x2 = x*x;
+ if (x2 > 1600) { // |x| > 40
+ double y2 = y*y;
+ if (x2 > 25e14) {// |x| > 5e7
+ double xy2 = (x*y)*(x*y);
+ return C((0.5 + y2 * (0.5 + 0.25*y2
+ - 0.16666666666666666667*xy2)) / x,
+ y * (-1 + y2 * (-0.66666666666666666667
+ + 0.13333333333333333333*xy2
+ - 0.26666666666666666667*y2))
+ / (2*x2 - 1));
+ }
+ return complex_mul_rc((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 D = spi2 * im_w_of_x(x);
+ double y2 = y*y;
+ return C
+ (D + y2 * (D + x - 2*D*x2)
+ + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2))
+ + x * (0.83333333333333333333
+ - 0.33333333333333333333 * x2)),
+ y * (1 - 2*D*x
+ + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2))
+ + y2*y2 * (0.26666666666666666667 -
+ x2 * (0.6 - 0.13333333333333333333 * x2)
+ - D*x * (1 - x2 * (1.3333333333333333333
+ - 0.26666666666666666667 * x2)))));
+ }
+ }
+} // cdawson
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-16 02:52:48 +0000