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