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, write to the Free Software
  53     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  54
  55
  56 #include "math.h"
  57 #include "math_private.h"
  58 #include <gnu/option-groups.h>
  59
  60 #if __OPTION_EGLIBC_LIBM_BIG
  61
  62 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
  63  * 1/sqrt(2) <= 1+x < sqrt(2)
  64  * Theoretical peak relative error = 5.3e-37,
  65  * relative peak error spread = 2.3e-14
  66  */
  67 static const long double
  68   P12 = 1.538612243596254322971797716843006400388E-6L,
  69   P11 = 4.998469661968096229986658302195402690910E-1L,
  70   P10 = 2.321125933898420063925789532045674660756E1L,
  71   P9 = 4.114517881637811823002128927449878962058E2L,
  72   P8 = 3.824952356185897735160588078446136783779E3L,
  73   P7 = 2.128857716871515081352991964243375186031E4L,
  74   P6 = 7.594356839258970405033155585486712125861E4L,
  75   P5 = 1.797628303815655343403735250238293741397E5L,
  76   P4 = 2.854829159639697837788887080758954924001E5L,
  77   P3 = 3.007007295140399532324943111654767187848E5L,
  78   P2 = 2.014652742082537582487669938141683759923E5L,
  79   P1 = 7.771154681358524243729929227226708890930E4L,
  80   P0 = 1.313572404063446165910279910527789794488E4L,
  81   /* Q12 = 1.000000000000000000000000000000000000000E0L, */
  82   Q11 = 4.839208193348159620282142911143429644326E1L,
  83   Q10 = 9.104928120962988414618126155557301584078E2L,
  84   Q9 = 9.147150349299596453976674231612674085381E3L,
  85   Q8 = 5.605842085972455027590989944010492125825E4L,
  86   Q7 = 2.248234257620569139969141618556349415120E5L,
  87   Q6 = 6.132189329546557743179177159925690841200E5L,
  88   Q5 = 1.158019977462989115839826904108208787040E6L,
  89   Q4 = 1.514882452993549494932585972882995548426E6L,
  90   Q3 = 1.347518538384329112529391120390701166528E6L,
  91   Q2 = 7.777690340007566932935753241556479363645E5L,
  92   Q1 = 2.626900195321832660448791748036714883242E5L,
  93   Q0 = 3.940717212190338497730839731583397586124E4L;
  94
  95 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
  96  * where z = 2(x-1)/(x+1)
  97  * 1/sqrt(2) <= x < sqrt(2)
  98  * Theoretical peak relative error = 1.1e-35,
  99  * relative peak error spread 1.1e-9
 100  */
 101 static const long double
 102   R5 = -8.828896441624934385266096344596648080902E-1L,
 103   R4 = 8.057002716646055371965756206836056074715E1L,
 104   R3 = -2.024301798136027039250415126250455056397E3L,
 105   R2 = 2.048819892795278657810231591630928516206E4L,
 106   R1 = -8.977257995689735303686582344659576526998E4L,
 107   R0 = 1.418134209872192732479751274970992665513E5L,
 108   /* S6 = 1.000000000000000000000000000000000000000E0L, */
 109   S5 = -1.186359407982897997337150403816839480438E2L,
 110   S4 = 3.998526750980007367835804959888064681098E3L,
 111   S3 = -5.748542087379434595104154610899551484314E4L,
 112   S2 = 4.001557694070773974936904547424676279307E5L,
 113   S1 = -1.332535117259762928288745111081235577029E6L,
 114   S0 = 1.701761051846631278975701529965589676574E6L;
 115
 116 /* C1 + C2 = ln 2 */
 117 static const long double C1 = 6.93145751953125E-1L;
 118 static const long double C2 = 1.428606820309417232121458176568075500134E-6L;
 119
 120 static const long double sqrth = 0.7071067811865475244008443621048490392848L;
 121 /* ln (2^16384 * (1 - 2^-113)) */
 122 static const long double maxlog = 1.1356523406294143949491931077970764891253E4L;
 123 static const long double zero = 0.0L;
 124
 125 long double
 126 __log1pl (long double xm1)
 127 {
 128   long double x, y, z, r, s;
 129   ieee854_long_double_shape_type u;
 130   int32_t hx;
 131   int e;
 132
 133   /* Test for NaN or infinity input. */
 134   u.value = xm1;
 135   hx = u.parts32.w0;
 136   if (hx >= 0x7fff0000)
 137     return xm1;
 138
 139   /* log1p(+- 0) = +- 0.  */
 140   if (((hx & 0x7fffffff) == 0)
 141       && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
 142     return xm1;
 143
 144   x = xm1 + 1.0L;
 145
 146   /* log1p(-1) = -inf */
 147   if (x <= 0.0L)
 148     {
 149       if (x == 0.0L)
 150         return (-1.0L / (x - x));
 151       else
 152         return (zero / (x - x));
 153     }
 154
 155   /* Separate mantissa from exponent.  */
 156
 157   /* Use frexp used so that denormal numbers will be handled properly.  */
 158   x = __frexpl (x, &e);
 159
 160   /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
 161      where z = 2(x-1)/x+1).  */
 162   if ((e > 2) || (e < -2))
 163     {
 164       if (x < sqrth)
 165         {                       /* 2( 2x-1 )/( 2x+1 ) */
 166           e -= 1;
 167           z = x - 0.5L;
 168           y = 0.5L * z + 0.5L;
 169         }
 170       else
 171         {                       /*  2 (x-1)/(x+1)   */
 172           z = x - 0.5L;
 173           z -= 0.5L;
 174           y = 0.5L * x + 0.5L;
 175         }
 176       x = z / y;
 177       z = x * x;
 178       r = ((((R5 * z
 179               + R4) * z
 180              + R3) * z
 181             + R2) * z
 182            + R1) * z
 183         + R0;
 184       s = (((((z
 185                + S5) * z
 186               + S4) * z
 187              + S3) * z
 188             + S2) * z
 189            + S1) * z
 190         + S0;
 191       z = x * (z * r / s);
 192       z = z + e * C2;
 193       z = z + x;
 194       z = z + e * C1;
 195       return (z);
 196     }
 197
 198
 199   /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
 200
 201   if (x < sqrth)
 202     {
 203       e -= 1;
 204       if (e != 0)
 205         x = 2.0L * x - 1.0L;    /*  2x - 1  */
 206       else
 207         x = xm1;
 208     }
 209   else
 210     {
 211       if (e != 0)
 212         x = x - 1.0L;
 213       else
 214         x = xm1;
 215     }
 216   z = x * x;
 217   r = (((((((((((P12 * x
 218                  + P11) * x
 219                 + P10) * x
 220                + P9) * x
 221               + P8) * x
 222              + P7) * x
 223             + P6) * x
 224            + P5) * x
 225           + P4) * x
 226          + P3) * x
 227         + P2) * x
 228        + P1) * x
 229     + P0;
 230   s = (((((((((((x
 231                  + Q11) * x
 232                 + Q10) * x
 233                + Q9) * x
 234               + Q8) * x
 235              + Q7) * x
 236             + Q6) * x
 237            + Q5) * x
 238           + Q4) * x
 239          + Q3) * x
 240         + Q2) * x
 241        + Q1) * x
 242     + Q0;
 243   y = x * (z * r / s);
 244   y = y + e * C2;
 245   z = y - 0.5L * z;
 246   z = z + x;
 247   z = z + e * C1;
 248   return (z);
 249 }
 250
 251 #else /* !__OPTION_EGLIBC_LIBM_BIG */
 252 # include <sysdeps/ieee754/ldbl-wrap/s_log1pl-wrap.c>
 253 #endif /* __OPTION_EGLIBC_LIBM_BIG */
 254
 255 weak_alias (__log1pl, log1pl)