winsup/cygwin/math/erfl.c

   1 /**
   2  * This file has no copyright assigned and is placed in the Public Domain.
   3  * This file is part of the mingw-w64 runtime package.
   4  * No warranty is given; refer to the file DISCLAIMER.PD within this package.
   5  */
   6 /*                                                      erfl.c
   7  *
   8  *      Error function
   9  *
  10  *
  11  *
  12  * SYNOPSIS:
  13  *
  14  * long double x, y, erfl();
  15  *
  16  * y = erfl( x );
  17  *
  18  *
  19  *
  20  * DESCRIPTION:
  21  *
  22  * The integral is
  23  *
  24  *                           x
  25  *                            -
  26  *                 2         | |          2
  27  *   erf(x)  =  --------     |    exp( - t  ) dt.
  28  *              sqrt(pi)   | |
  29  *                          -
  30  *                           0
  31  *
  32  * The magnitude of x is limited to about 106.56 for IEEE
  33  * arithmetic; 1 or -1 is returned outside this range.
  34  *
  35  * For 0 <= |x| < 1, erf(x) = x * P6(x^2)/Q6(x^2);
  36  * Otherwise: erf(x) = 1 - erfc(x).
  37  *
  38  *
  39  *
  40  * ACCURACY:
  41  *
  42  *                      Relative error:
  43  * arithmetic   domain     # trials      peak         rms
  44  *    IEEE      0,1         50000       2.0e-19     5.7e-20
  45  *
  46  */
  47
  48 /*                                                      erfcl.c
  49  *
  50  *      Complementary error function
  51  *
  52  *
  53  *
  54  * SYNOPSIS:
  55  *
  56  * long double x, y, erfcl();
  57  *
  58  * y = erfcl( x );
  59  *
  60  *
  61  *
  62  * DESCRIPTION:
  63  *
  64  *
  65  *  1 - erf(x) =
  66  *
  67  *                           inf.
  68  *                             -
  69  *                  2         | |          2
  70  *   erfc(x)  =  --------     |    exp( - t  ) dt
  71  *               sqrt(pi)   | |
  72  *                           -
  73  *                            x
  74  *
  75  *
  76  * For small x, erfc(x) = 1 - erf(x); otherwise rational
  77  * approximations are computed.
  78  *
  79  * A special function expx2l.c is used to suppress error amplification
  80  * in computing exp(-x^2).
  81  *
  82  *
  83  * ACCURACY:
  84  *
  85  *                      Relative error:
  86  * arithmetic   domain     # trials      peak         rms
  87  *    IEEE      0,13        50000      8.4e-19      9.7e-20
  88  *    IEEE      6,106.56    20000      2.9e-19      7.1e-20
  89  *
  90  *
  91  * ERROR MESSAGES:
  92  *
  93  *   message          condition              value returned
  94  * erfcl underflow    x^2 > MAXLOGL              0.0
  95  *
  96  *
  97  */
  98
  99
 100 /*
 101 Modified from file ndtrl.c
 102 Cephes Math Library Release 2.3:  January, 1995
 103 Copyright 1984, 1995 by Stephen L. Moshier
 104 */
 105
 106 #include <math.h>
 107 #include "cephes_mconf.h"
 108
 109 long double erfl(long double x);
 110
 111 /* erfc(x) = exp(-x^2) P(1/x)/Q(1/x)
 112    1/8 <= 1/x <= 1
 113    Peak relative error 5.8e-21  */
 114
 115 static const uLD P[10] = {
 116   { { 0x4bf0,0x9ad8,0x7a03,0x86c7,0x401d, 0, 0, 0 } },
 117   { { 0xdf23,0xd843,0x4032,0x8881,0x401e, 0, 0, 0 } },
 118   { { 0xd025,0xcfd5,0x8494,0x88d3,0x401e, 0, 0, 0 } },
 119   { { 0xb6d0,0xc92b,0x5417,0xacb1,0x401d, 0, 0, 0 } },
 120   { { 0xada8,0x356a,0x4982,0x94a6,0x401c, 0, 0, 0 } },
 121   { { 0x4e13,0xcaee,0x9e31,0xb258,0x401a, 0, 0, 0 } },
 122   { { 0x5840,0x554d,0x37a3,0x9239,0x4018, 0, 0, 0 } },
 123   { { 0x3b58,0x3da2,0xaf02,0x9780,0x4015, 0, 0, 0 } },
 124   { { 0x0144,0x489e,0xbe68,0x9c31,0x4011, 0, 0, 0 } },
 125   { { 0x333b,0xd9e6,0xd404,0x986f,0xbfee, 0, 0, 0 } }
 126 };
 127 static const uLD Q[] = {
 128   { { 0x0e43,0x302d,0x79ed,0x86c7,0x401d, 0, 0, 0 } },
 129   { { 0xf817,0x9128,0xc0f8,0xd48b,0x401e, 0, 0, 0 } },
 130   { { 0x8eae,0x8dad,0x6eb4,0x9aa2,0x401f, 0, 0, 0 } },
 131   { { 0x00e7,0x7595,0xcd06,0x88bb,0x401f, 0, 0, 0 } },
 132   { { 0x4991,0xcfda,0x52f1,0xa2a9,0x401e, 0, 0, 0 } },
 133   { { 0xc39d,0xe415,0xc43d,0x87c0,0x401d, 0, 0, 0 } },
 134   { { 0xa75d,0x436f,0x30dd,0xa027,0x401b, 0, 0, 0 } },
 135   { { 0xc4cb,0x305a,0xbf78,0x8220,0x4019, 0, 0, 0 } },
 136   { { 0x3708,0x33b1,0x07fa,0x8644,0x4016, 0, 0, 0 } },
 137   { { 0x24fa,0x96f6,0x7153,0x8a6c,0x4012, 0, 0, 0 } }
 138 };
 139
 140 /* erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
 141    1/128 <= 1/x < 1/8
 142    Peak relative error 1.9e-21  */
 143
 144 static const uLD R[] = {
 145   { { 0x260a,0xab95,0x2fc7,0xe7c4,0x4000, 0, 0, 0 } },
 146   { { 0x4761,0x613e,0xdf6d,0xe58e,0x4001, 0, 0, 0 } },
 147   { { 0x0615,0x4b00,0x575f,0xdc7b,0x4000, 0, 0, 0 } },
 148   { { 0x521d,0x8527,0x3435,0x8dc2,0x3ffe, 0, 0, 0 } },
 149   { { 0x22cf,0xc711,0x6c5b,0xdcfb,0x3ff9, 0, 0, 0 } }
 150 };
 151 static const uLD S[] = {
 152   { { 0x5de6,0x17d7,0x54d6,0xaba9,0x4002, 0, 0, 0 } },
 153   { { 0x55d5,0xd300,0xe71e,0xf564,0x4002, 0, 0, 0 } },
 154   { { 0xb611,0x8f76,0xf020,0xd255,0x4001, 0, 0, 0 } },
 155   { { 0x3684,0x3798,0xb793,0x80b0,0x3fff, 0, 0, 0 } },
 156   { { 0xf5af,0x2fb2,0x1e57,0xc3d7,0x3ffa, 0, 0, 0 } }
 157 };
 158
 159 /* erf(x)  = x T(x^2)/U(x^2)
 160    0 <= x <= 1
 161    Peak relative error 7.6e-23  */
 162
 163 static const uLD T[] = {
 164   { { 0xfd7a,0x3a1a,0x705b,0xe0c4,0x3ffb, 0, 0, 0 } },
 165   { { 0x3128,0xc337,0x3716,0xace5,0x4001, 0, 0, 0 } },
 166   { { 0x9517,0x4e93,0x540e,0x8f97,0x4007, 0, 0, 0 } },
 167   { { 0x6118,0x6059,0x9093,0xa757,0x400a, 0, 0, 0 } },
 168   { { 0xb954,0xa987,0xc60c,0xbc83,0x400e, 0, 0, 0 } },
 169   { { 0x7a56,0xe45a,0xa4bd,0x975b,0x4010, 0, 0, 0 } },
 170   { { 0xc446,0x6bab,0x0b2a,0x86d0,0x4013, 0, 0, 0 } }
 171 };
 172
 173 static const uLD U[] = {
 174   { { 0x3453,0x1f8e,0xf688,0xb507,0x4004, 0, 0, 0 } },
 175   { { 0x71ac,0xb12f,0x21ca,0xf2e2,0x4008, 0, 0, 0 } },
 176   { { 0xffe8,0x9cac,0x3b84,0xc2ac,0x400c, 0, 0, 0 } },
 177   { { 0x481d,0x445b,0xc807,0xc232,0x400f, 0, 0, 0 } },
 178   { { 0x9ad5,0x1aef,0x45b1,0xe25e,0x4011, 0, 0, 0 } },
 179   { { 0x71a7,0x1cad,0x012e,0xeef3,0x4012, 0, 0, 0 } }
 180 };
 181
 182 /*                                                      expx2l.c
 183  *
 184  *      Exponential of squared argument
 185  *
 186  *
 187  *
 188  * SYNOPSIS:
 189  *
 190  * long double x, y, expmx2l();
 191  * int sign;
 192  *
 193  * y = expx2l( x );
 194  *
 195  *
 196  *
 197  * DESCRIPTION:
 198  *
 199  * Computes y = exp(x*x) while suppressing error amplification
 200  * that would ordinarily arise from the inexactness of the
 201  * exponential argument x*x.
 202  *
 203  *
 204  *
 205  * ACCURACY:
 206  *
 207  *                      Relative error:
 208  * arithmetic      domain        # trials      peak         rms
 209  *   IEEE     -106.566, 106.566    10^5       1.6e-19     4.4e-20
 210  *
 211  */
 212
 213 #define M 32768.0L
 214 #define MINV 3.0517578125e-5L
 215
 216 static long double expx2l (long double x)
 217 {
 218         long double u, u1, m, f;
 219
 220         x = fabsl (x);
 221         /* Represent x as an exact multiple of M plus a residual.
 222            M is a power of 2 chosen so that exp(m * m) does not overflow
 223            or underflow and so that |x - m| is small.  */
 224         m = MINV * floorl(M * x + 0.5L);
 225         f = x - m;
 226
 227         /* x^2 = m^2 + 2mf + f^2 */
 228         u = m * m;
 229         u1 = 2 * m * f  +  f * f;
 230
 231         if ((u + u1) > MAXLOGL)
 232                 return (INFINITYL);
 233
 234         /* u is exact, u1 is small.  */
 235         u = expl(u) * expl(u1);
 236         return (u);
 237 }
 238
 239 long double erfcl(long double a)
 240 {
 241         long double p, q, x, y, z;
 242
 243         if (isinf (a))
 244                 return (signbit(a) ? 2.0 : 0.0);
 245
 246         if (isnan (a))
 247                 return (a);
 248
 249         x = fabsl (a);
 250
 251         if (x < 1.0L)
 252                 return (1.0L - erfl(a));
 253
 254         z = a * a;
 255
 256         if (z  > MAXLOGL)
 257         {
 258 under:
 259                 mtherr("erfcl", UNDERFLOW);
 260                 errno = ERANGE;
 261                 return (signbit(a) ? 2.0 : 0.0);
 262         }
 263
 264         /* Compute z = expl(a * a).  */
 265         z = expx2l(a);
 266         y = 1.0L/x;
 267
 268         if (x < 8.0L)
 269         {
 270                 p = polevll(y, P, 9);
 271                 q = p1evll(y, Q, 10);
 272         }
 273         else
 274         {
 275                 q = y * y;
 276                 p = y * polevll(q, R, 4);
 277                 q = p1evll(q, S, 5);
 278         }
 279         y = p/(q * z);
 280
 281         if (a < 0.0L)
 282                 y = 2.0L - y;
 283
 284         if (y == 0.0L)
 285                 goto under;
 286
 287         return (y);
 288 }
 289
 290 long double erfl(long double x)
 291 {
 292         long double y, z;
 293
 294         if (x == 0.0L)
 295                 return (x);
 296
 297         if (isinf (x))
 298                 return (signbit(x) ?  -1.0L : 1.0L);
 299
 300         if (fabsl(x) > 1.0L)
 301                 return (1.0L - erfcl(x));
 302
 303         z = x * x;
 304         y = x * polevll(z, T, 6) / p1evll(z, U, 6);
 305         return (y);
 306 }