sysdeps/ieee754/ldbl-128/s_log1pl.c

   1 /*                                                      log1pl.c
   2  *
   3  *      Relative error logarithm
   4  *      Natural logarithm of 1+x, 128-bit long double precision
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * long double x, y, log1pl();
  11  *
  12  * y = log1pl( x );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns the base e (2.718...) logarithm of 1+x.
  19  *
  20  * The argument 1+x is separated into its exponent and fractional
  21  * parts.  If the exponent is between -1 and +1, the logarithm
  22  * of the fraction is approximated by
  23  *
  24  *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
  25  *
  26  * Otherwise, setting  z = 2(w-1)/(w+1),
  27  *
  28  *     log(w) = z + z^3 P(z)/Q(z).
  29  *
  30  *
  31  *
  32  * ACCURACY:
  33  *
  34  *                      Relative error:
  35  * arithmetic   domain     # trials      peak         rms
  36  *    IEEE      -1, 8       100000      1.9e-34     4.3e-35
  37  */
  38
  39 /* Copyright 2001 by Stephen L. Moshier
  40
  41     This library is free software; you can redistribute it and/or
  42     modify it under the terms of the GNU Lesser General Public
  43     License as published by the Free Software Foundation; either
  44     version 2.1 of the License, or (at your option) any later version.
  45
  46     This library is distributed in the hope that it will be useful,
  47     but WITHOUT ANY WARRANTY; without even the implied warranty of
  48     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  49     Lesser General Public License for more details.
  50
  51     You should have received a copy of the GNU Lesser General Public
  52     License along with this library; if not, see
  53     <http://www.gnu.org/licenses/>.  */
  54
  55
  56 #include <math.h>
  57 #include <math_private.h>
  58
  59 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
  60  * 1/sqrt(2) <= 1+x < sqrt(2)
  61  * Theoretical peak relative error = 5.3e-37,
  62  * relative peak error spread = 2.3e-14
  63  */
  64 static const long double
  65   P12 = 1.538612243596254322971797716843006400388E-6L,
  66   P11 = 4.998469661968096229986658302195402690910E-1L,
  67   P10 = 2.321125933898420063925789532045674660756E1L,
  68   P9 = 4.114517881637811823002128927449878962058E2L,
  69   P8 = 3.824952356185897735160588078446136783779E3L,
  70   P7 = 2.128857716871515081352991964243375186031E4L,
  71   P6 = 7.594356839258970405033155585486712125861E4L,
  72   P5 = 1.797628303815655343403735250238293741397E5L,
  73   P4 = 2.854829159639697837788887080758954924001E5L,
  74   P3 = 3.007007295140399532324943111654767187848E5L,
  75   P2 = 2.014652742082537582487669938141683759923E5L,
  76   P1 = 7.771154681358524243729929227226708890930E4L,
  77   P0 = 1.313572404063446165910279910527789794488E4L,
  78   /* Q12 = 1.000000000000000000000000000000000000000E0L, */
  79   Q11 = 4.839208193348159620282142911143429644326E1L,
  80   Q10 = 9.104928120962988414618126155557301584078E2L,
  81   Q9 = 9.147150349299596453976674231612674085381E3L,
  82   Q8 = 5.605842085972455027590989944010492125825E4L,
  83   Q7 = 2.248234257620569139969141618556349415120E5L,
  84   Q6 = 6.132189329546557743179177159925690841200E5L,
  85   Q5 = 1.158019977462989115839826904108208787040E6L,
  86   Q4 = 1.514882452993549494932585972882995548426E6L,
  87   Q3 = 1.347518538384329112529391120390701166528E6L,
  88   Q2 = 7.777690340007566932935753241556479363645E5L,
  89   Q1 = 2.626900195321832660448791748036714883242E5L,
  90   Q0 = 3.940717212190338497730839731583397586124E4L;
  91
  92 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
  93  * where z = 2(x-1)/(x+1)
  94  * 1/sqrt(2) <= x < sqrt(2)
  95  * Theoretical peak relative error = 1.1e-35,
  96  * relative peak error spread 1.1e-9
  97  */
  98 static const long double
  99   R5 = -8.828896441624934385266096344596648080902E-1L,
 100   R4 = 8.057002716646055371965756206836056074715E1L,
 101   R3 = -2.024301798136027039250415126250455056397E3L,
 102   R2 = 2.048819892795278657810231591630928516206E4L,
 103   R1 = -8.977257995689735303686582344659576526998E4L,
 104   R0 = 1.418134209872192732479751274970992665513E5L,
 105   /* S6 = 1.000000000000000000000000000000000000000E0L, */
 106   S5 = -1.186359407982897997337150403816839480438E2L,
 107   S4 = 3.998526750980007367835804959888064681098E3L,
 108   S3 = -5.748542087379434595104154610899551484314E4L,
 109   S2 = 4.001557694070773974936904547424676279307E5L,
 110   S1 = -1.332535117259762928288745111081235577029E6L,
 111   S0 = 1.701761051846631278975701529965589676574E6L;
 112
 113 /* C1 + C2 = ln 2 */
 114 static const long double C1 = 6.93145751953125E-1L;
 115 static const long double C2 = 1.428606820309417232121458176568075500134E-6L;
 116
 117 static const long double sqrth = 0.7071067811865475244008443621048490392848L;
 118 /* ln (2^16384 * (1 - 2^-113)) */
 119 static const long double maxlog = 1.1356523406294143949491931077970764891253E4L;
 120 static const long double zero = 0.0L;
 121
 122 long double
 123 __log1pl (long double xm1)
 124 {
 125   long double x, y, z, r, s;
 126   ieee854_long_double_shape_type u;
 127   int32_t hx;
 128   int e;
 129
 130   /* Test for NaN or infinity input. */
 131   u.value = xm1;
 132   hx = u.parts32.w0;
 133   if (hx >= 0x7fff0000)
 134     return xm1;
 135
 136   /* log1p(+- 0) = +- 0.  */
 137   if (((hx & 0x7fffffff) == 0)
 138       && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
 139     return xm1;
 140
 141   if ((hx & 0x7fffffff) < 0x3f8e0000)
 142     {
 143       if ((int) xm1 == 0)
 144         return xm1;
 145     }
 146
 147   x = xm1 + 1.0L;
 148
 149   /* log1p(-1) = -inf */
 150   if (x <= 0.0L)
 151     {
 152       if (x == 0.0L)
 153         return (-1.0L / (x - x));
 154       else
 155         return (zero / (x - x));
 156     }
 157
 158   /* Separate mantissa from exponent.  */
 159
 160   /* Use frexp used so that denormal numbers will be handled properly.  */
 161   x = __frexpl (x, &e);
 162
 163   /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
 164      where z = 2(x-1)/x+1).  */
 165   if ((e > 2) || (e < -2))
 166     {
 167       if (x < sqrth)
 168         {                       /* 2( 2x-1 )/( 2x+1 ) */
 169           e -= 1;
 170           z = x - 0.5L;
 171           y = 0.5L * z + 0.5L;
 172         }
 173       else
 174         {                       /*  2 (x-1)/(x+1)   */
 175           z = x - 0.5L;
 176           z -= 0.5L;
 177           y = 0.5L * x + 0.5L;
 178         }
 179       x = z / y;
 180       z = x * x;
 181       r = ((((R5 * z
 182               + R4) * z
 183              + R3) * z
 184             + R2) * z
 185            + R1) * z
 186         + R0;
 187       s = (((((z
 188                + S5) * z
 189               + S4) * z
 190              + S3) * z
 191             + S2) * z
 192            + S1) * z
 193         + S0;
 194       z = x * (z * r / s);
 195       z = z + e * C2;
 196       z = z + x;
 197       z = z + e * C1;
 198       return (z);
 199     }
 200
 201
 202   /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
 203
 204   if (x < sqrth)
 205     {
 206       e -= 1;
 207       if (e != 0)
 208         x = 2.0L * x - 1.0L;    /*  2x - 1  */
 209       else
 210         x = xm1;
 211     }
 212   else
 213     {
 214       if (e != 0)
 215         x = x - 1.0L;
 216       else
 217         x = xm1;
 218     }
 219   z = x * x;
 220   r = (((((((((((P12 * x
 221                  + P11) * x
 222                 + P10) * x
 223                + P9) * x
 224               + P8) * x
 225              + P7) * x
 226             + P6) * x
 227            + P5) * x
 228           + P4) * x
 229          + P3) * x
 230         + P2) * x
 231        + P1) * x
 232     + P0;
 233   s = (((((((((((x
 234                  + Q11) * x
 235                 + Q10) * x
 236                + Q9) * x
 237               + Q8) * x
 238              + Q7) * x
 239             + Q6) * x
 240            + Q5) * x
 241           + Q4) * x
 242          + Q3) * x
 243         + Q2) * x
 244        + Q1) * x
 245     + Q0;
 246   y = x * (z * r / s);
 247   y = y + e * C2;
 248   z = y - 0.5L * z;
 249   z = z + x;
 250   z = z + e * C1;
 251   return (z);
 252 }
 253
 254 weak_alias (__log1pl, log1pl)