libquadmath/math/log1pq.c

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