From 02626f616d01a8166f5a7c275a2d1fd5aeb150f9 Mon Sep 17 00:00:00 2001 From: Danny Smith Date: Sun, 8 May 2005 08:07:18 +0000 Subject: [PATCH] * mingwex/math/erfl.c: New file. * mingwex/Makefile.i (MATH_DISTFILES): Add erfl.c. (MATH_OBJS): Add erfl.o. * include/math.h (erfl, erfcl): Uncomment prototypes. --- winsup/mingw/ChangeLog | 7 + winsup/mingw/include/math.h | 4 - winsup/mingw/mingwex/Makefile.in | 4 +- winsup/mingw/mingwex/math/erfl.c | 299 +++++++++++++++++++++++++++++++ 4 files changed, 308 insertions(+), 6 deletions(-) create mode 100755 winsup/mingw/mingwex/math/erfl.c diff --git a/winsup/mingw/ChangeLog b/winsup/mingw/ChangeLog index 989da3b19..ac21c8c99 100644 --- a/winsup/mingw/ChangeLog +++ b/winsup/mingw/ChangeLog @@ -1,3 +1,10 @@ +2005-05-08 Danny Smith + + * mingwex/math/erfl.c: New file. + * mingwex/Makefile.i (MATH_DISTFILES): Add erfl.c. + (MATH_OBJS): Add erfl.o. + * include/math.h (erfl, erfcl): Uncomment prototypes. + 2005-05-04 Danny Smith * include/wchar.h (WCHAR_MAX): Define as 0xffff, so preprocessor diff --git a/winsup/mingw/include/math.h b/winsup/mingw/include/math.h index 585723b22..751fa9509 100644 --- a/winsup/mingw/include/math.h +++ b/winsup/mingw/include/math.h @@ -565,16 +565,12 @@ extern long double __cdecl sqrtl (long double); /* 7.12.8.1 The erf functions */ extern double __cdecl erf (double); extern float __cdecl erff (float); -/* TODO extern long double __cdecl erfl (long double); -*/ /* 7.12.8.2 The erfc functions */ extern double __cdecl erfc (double); extern float __cdecl erfcf (float); -/* TODO extern long double __cdecl erfcl (long double); -*/ /* 7.12.8.3 The lgamma functions */ extern double __cdecl lgamma (double); diff --git a/winsup/mingw/mingwex/Makefile.in b/winsup/mingw/mingwex/Makefile.in index 8876cf9bd..8bfebc533 100644 --- a/winsup/mingw/mingwex/Makefile.in +++ b/winsup/mingw/mingwex/Makefile.in @@ -41,7 +41,7 @@ MATH_DISTFILES = \ atanf.c atanl.c cbrt.c cbrtf.c cbrtl.c ceilf.S ceill.S \ cephes_emath.h cephes_emath.c cephes_mconf.h \ copysign.S copysignf.S copysignl.S cosf.S coshf.c coshl.c cosl.S \ - exp2.S exp2f.S exp2l.S expf.c expl.c expm1.c expm1l.c expm1f.c \ + erfl.c exp2.S exp2f.S exp2l.S expf.c expl.c expm1.c expm1l.c expm1f.c \ fabs.c fabsf.c fabsl.c \ fdim.c fdimf.c fdiml.c floorf.S floorl.S fma.S fmaf.S fmal.c \ fmax.c fmaxf.c fmaxl.c fmin.c fminf.c fminl.c fmodf.c \ @@ -133,7 +133,7 @@ MATH_OBJS = \ atanf.o atanl.o cbrt.o cbrtf.o cbrtl.o ceilf.o ceill.o \ cephes_emath.o \ copysign.o copysignf.o copysignl.o cosf.o coshf.o coshl.o cosl.o \ - exp2.o exp2f.o exp2l.o expf.o expl.o expm1.o expm1l.o expm1f.o \ + erfl.o exp2.o exp2f.o exp2l.o expf.o expl.o expm1.o expm1l.o expm1f.o \ fabs.o fabsf.o fabsl.o \ fdim.o fdimf.o fdiml.o floorf.o floorl.o fma.o fmaf.o fmal.o \ fmax.o fmaxf.o fmaxl.o fmin.o fminf.o fminl.o fmodf.o \ diff --git a/winsup/mingw/mingwex/math/erfl.c b/winsup/mingw/mingwex/math/erfl.c new file mode 100755 index 000000000..ac50d98c2 --- /dev/null +++ b/winsup/mingw/mingwex/math/erfl.c @@ -0,0 +1,299 @@ +/* erfl.c + * + * Error function + * + * + * + * SYNOPSIS: + * + * long double x, y, erfl(); + * + * y = erfl( x ); + * + * + * + * DESCRIPTION: + * + * The integral is + * + * x + * - + * 2 | | 2 + * erf(x) = -------- | exp( - t ) dt. + * sqrt(pi) | | + * - + * 0 + * + * The magnitude of x is limited to about 106.56 for IEEE + * arithmetic; 1 or -1 is returned outside this range. + * + * For 0 <= |x| < 1, erf(x) = x * P6(x^2)/Q6(x^2); + * Otherwise: erf(x) = 1 - erfc(x). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,1 50000 2.0e-19 5.7e-20 + * + */ + +/* erfcl.c + * + * Complementary error function + * + * + * + * SYNOPSIS: + * + * long double x, y, erfcl(); + * + * y = erfcl( x ); + * + * + * + * DESCRIPTION: + * + * + * 1 - erf(x) = + * + * inf. + * - + * 2 | | 2 + * erfc(x) = -------- | exp( - t ) dt + * sqrt(pi) | | + * - + * x + * + * + * For small x, erfc(x) = 1 - erf(x); otherwise rational + * approximations are computed. + * + * A special function expx2l.c is used to suppress error amplification + * in computing exp(-x^2). + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,13 50000 8.4e-19 9.7e-20 + * IEEE 6,106.56 20000 2.9e-19 7.1e-20 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * erfcl underflow x^2 > MAXLOGL 0.0 + * + * + */ + + +/* +Modified from file ndtrl.c +Cephes Math Library Release 2.3: January, 1995 +Copyright 1984, 1995 by Stephen L. Moshier +*/ + +#include +#include "cephes_mconf.h" + +/* erfc(x) = exp(-x^2) P(1/x)/Q(1/x) + 1/8 <= 1/x <= 1 + Peak relative error 5.8e-21 */ + +static const unsigned short P[] = { +0x4bf0,0x9ad8,0x7a03,0x86c7,0x401d, XPD +0xdf23,0xd843,0x4032,0x8881,0x401e, XPD +0xd025,0xcfd5,0x8494,0x88d3,0x401e, XPD +0xb6d0,0xc92b,0x5417,0xacb1,0x401d, XPD +0xada8,0x356a,0x4982,0x94a6,0x401c, XPD +0x4e13,0xcaee,0x9e31,0xb258,0x401a, XPD +0x5840,0x554d,0x37a3,0x9239,0x4018, XPD +0x3b58,0x3da2,0xaf02,0x9780,0x4015, XPD +0x0144,0x489e,0xbe68,0x9c31,0x4011, XPD +0x333b,0xd9e6,0xd404,0x986f,0xbfee, XPD +}; +static const unsigned short Q[] = { +/* 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD */ +0x0e43,0x302d,0x79ed,0x86c7,0x401d, XPD +0xf817,0x9128,0xc0f8,0xd48b,0x401e, XPD +0x8eae,0x8dad,0x6eb4,0x9aa2,0x401f, XPD +0x00e7,0x7595,0xcd06,0x88bb,0x401f, XPD +0x4991,0xcfda,0x52f1,0xa2a9,0x401e, XPD +0xc39d,0xe415,0xc43d,0x87c0,0x401d, XPD +0xa75d,0x436f,0x30dd,0xa027,0x401b, XPD +0xc4cb,0x305a,0xbf78,0x8220,0x4019, XPD +0x3708,0x33b1,0x07fa,0x8644,0x4016, XPD +0x24fa,0x96f6,0x7153,0x8a6c,0x4012, XPD +}; + +/* erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2) + 1/128 <= 1/x < 1/8 + Peak relative error 1.9e-21 */ + +static const unsigned short R[] = { +0x260a,0xab95,0x2fc7,0xe7c4,0x4000, XPD +0x4761,0x613e,0xdf6d,0xe58e,0x4001, XPD +0x0615,0x4b00,0x575f,0xdc7b,0x4000, XPD +0x521d,0x8527,0x3435,0x8dc2,0x3ffe, XPD +0x22cf,0xc711,0x6c5b,0xdcfb,0x3ff9, XPD +}; +static const unsigned short S[] = { +/* 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD */ +0x5de6,0x17d7,0x54d6,0xaba9,0x4002, XPD +0x55d5,0xd300,0xe71e,0xf564,0x4002, XPD +0xb611,0x8f76,0xf020,0xd255,0x4001, XPD +0x3684,0x3798,0xb793,0x80b0,0x3fff, XPD +0xf5af,0x2fb2,0x1e57,0xc3d7,0x3ffa, XPD +}; + +/* erf(x) = x T(x^2)/U(x^2) + 0 <= x <= 1 + Peak relative error 7.6e-23 */ + +static const unsigned short T[] = { +0xfd7a,0x3a1a,0x705b,0xe0c4,0x3ffb, XPD +0x3128,0xc337,0x3716,0xace5,0x4001, XPD +0x9517,0x4e93,0x540e,0x8f97,0x4007, XPD +0x6118,0x6059,0x9093,0xa757,0x400a, XPD +0xb954,0xa987,0xc60c,0xbc83,0x400e, XPD +0x7a56,0xe45a,0xa4bd,0x975b,0x4010, XPD +0xc446,0x6bab,0x0b2a,0x86d0,0x4013, XPD +}; + +static const unsigned short U[] = { +/* 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD */ +0x3453,0x1f8e,0xf688,0xb507,0x4004, XPD +0x71ac,0xb12f,0x21ca,0xf2e2,0x4008, XPD +0xffe8,0x9cac,0x3b84,0xc2ac,0x400c, XPD +0x481d,0x445b,0xc807,0xc232,0x400f, XPD +0x9ad5,0x1aef,0x45b1,0xe25e,0x4011, XPD +0x71a7,0x1cad,0x012e,0xeef3,0x4012, XPD +}; + +/* expx2l.c + * + * Exponential of squared argument + * + * + * + * SYNOPSIS: + * + * long double x, y, expmx2l(); + * int sign; + * + * y = expx2l( x ); + * + * + * + * DESCRIPTION: + * + * Computes y = exp(x*x) while suppressing error amplification + * that would ordinarily arise from the inexactness of the + * exponential argument x*x. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -106.566, 106.566 10^5 1.6e-19 4.4e-20 + * + */ + +#define M 32768.0L +#define MINV 3.0517578125e-5L + +static long double expx2l (long double x) +{ + long double u, u1, m, f; + + x = fabsl (x); + /* Represent x as an exact multiple of M plus a residual. + M is a power of 2 chosen so that exp(m * m) does not overflow + or underflow and so that |x - m| is small. */ + m = MINV * floorl(M * x + 0.5L); + f = x - m; + + /* x^2 = m^2 + 2mf + f^2 */ + u = m * m; + u1 = 2 * m * f + f * f; + + if ((u+u1) > MAXLOGL) + return (INFINITYL); + + /* u is exact, u1 is small. */ + u = expl(u) * expl(u1); + return(u); +} + +long double erfcl(long double a) +{ +long double p,q,x,y,z; + +if (isinf (a)) + return (signbit (a) ? 2.0 : 0.0); + +x = fabsl (a); + +if (x < 1.0L) + return (1.0L - erfl(a)); + +z = a * a; + +if( z > MAXLOGL ) + { +under: + mtherr( "erfcl", UNDERFLOW ); + errno = ERANGE; + return (signbit (a) ? 2.0 : 0.0); + } + +/* Compute z = expl(a * a). */ +z = expx2l (a); +y = 1.0L/x; + +if (x < 8.0L) + { + p = polevll (y, P, 9); + q = p1evll (y, Q, 10); + } +else + { + q = y * y; + p = y * polevll (q, R, 4); + q = p1evll (q, S, 5); + } +y = p/(q * z); + +if (a < 0.0L) + y = 2.0L - y; + +if (y == 0.0L) + goto under; + +return (y); +} + +long double erfl(long double x) +{ +long double y, z; + +if( x == 0.0L ) + return (x); + +if (isinf (x)) + return (signbit (x) ? -1.0L : 1.0L); + +if (fabsl(x) > 1.0L) + return (1.0L - erfcl (x)); + +z = x * x; +y = x * polevll( z, T, 6 ) / p1evll( z, U, 6 ); +return( y ); +}